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Abstract 


Covariance  Sequence  Approximation  with  Applications 
to  Spectrun  Analysis  and  Digital  Filter  Design 


L 


Modeling  and  estimation  procedures  for  covariance  sequence  and 
spectrum  approximation  are  developed  in  this  thesis.  The  covariance 
sequence  is  modeled  as  a complex  linear  combination  of  damped  complex 
exponentials.  This  model  arises  naturally  as  the  representation  for 
the  covariance  sequence  associated  with  a strictly  proper  ARMA(M,N) 
system  driven  by  white  noise.  Related  to  this  seemingly  natural  covar- 
iance model  is  a synthesis  procedure  for  a subclass  of  wide  sense 
stationary  ARMA(M,N)  processes.  The  resulting  spectral  representation 
for  the  covariance  sequence  is  a positive  real  linear  combination  of 
damped  complex  exponentials,  a generalization  of  the  standard  represen- 
tation in  terms  of  stochastic  almost  periodic  functions.  The  importance 
of  the  generalized  structure  lies  in  the  more  efficient  representation 
of  a large  class  of  wide  sense  stationary  processes.  For  this 
parametric  approach  estimation  techniques  are  developed  that  lie  in 
philosophy  somewhere  between  the  nonlinear  least  squares  approach  and 
the  tractable  modified  least  squares  procedure.  The  resulting  parameter 
estimation  equations  are  linear,  except  for  a single  polynomial  root- 
f India g problem  that  must  be  solved. 

Parametric  spectra  are  efficiently  approximated  for  a large  class 
of  processes  using  the  covariance  sequence  approximant.  The  generality 
of  the  approximant  does  not  force  the  sequence  to  be  non-negative 
definite,  so  that  the  corresponding  "spectrum"  estimate  can  be  negative 
at  some  frequencies.  The  same  generality  of  the  model  permits  an 
accurate  representation  of  sequences  that  are  not  non-negative  definite. 
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Recursive  digital  filter  design  based  on  second  order  statistical 
information  is  achieved  by  covariance  sequence  approximation.  The 
stochastic  minimal  realization  problem  encountered  in  covariance 
invariant  recursive  digital  filter  design  forms  a special  case , and 
results  in  stable  ARMA.(M,M-1)  filters  with  real  coefficients.  For  a 
general  desired  spectrum  the  resulting  filter  is  stable  and  of  ARMA(2M, 
2tt-2)  type  with  real  coefficients.  The  latter  design  is  efficiently 
achieved  by  approximating  the  covariance  sequence  associated  with  the 
square  root  of  the  desired  spectrum.  The  phase  characteristic  associated 
with  these  designs  can  be  manipulated  to  some  extent  to  achieve  linear 
phase,  phase  equalization,  or  filter  lnvertibility  properties. 


1 INTRODUCTION 

The  theory  of  parametric  statistical  inference  is  concerned  with 
the  identification  of  parametric  models  from  random  data.  When  the  data 
arises  in  the  form  of  a time  series,  then  the  basic  problem  is  to  identi- 
fy a parametric  model  for  the  underlying  stochastic  process.  If  the 
underlying  process  is  normal,  then  the  covariance  sequence  provides  a 
complete  statistical  characterization  of  the  process.  Furthermore,  the 
covariance  sequence  for  a process  lies  in  one-to-one  correspondence  with 
the  power  spectral  density  and  summarizes  all  the  information  required 
to  realize  (or  approximate)  a minimum  mean-squared  error  (MMSE)  filter 
for  the  process.  If  the  power  spectral  density  happens  to  be  rational, 
then  the  MMSE  filter  is  of  Wiener/Kaltcan  variety. 

All  of  this  suggests  that  one  might  profitably  describe  a time 
series  by  identifying  a parametric  model  for  its  covariance  sequence. 

In  this  thesis  we  develop  techniques  for  identifying  parameters  in  the 
following  parametric  model  of  the  covariance  sequence: 


*k  = ; Ai  v 1 l-1 

This  model  is  a useful  one  because  it  includes  the  covariance  sequences 
of  the  following  processes  as  special  cases:  (i)  autoregressive  moving 
average  (ARMA) , (ii)  random  amplitude,  random  phase  sinusoid,  (iii)  white 
noise,  and  (iv)  linear  combinations  of  (i)-(iil).  Furthermore,  the 


model  reduces  under  special  circumstances  to  the  model  implicit  in  DFT 
analysis  of  random  data  and  in  Pisarenko's  decomposition  of  a covariance 


l 


I 


sequence. 
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Once  the  parameters  {A^},  {p^}  of  (1.1)  have  been  determined,  it  is 
a simple  matter  to  find  the  corresponding  parametric  spectral  density 
estimate.  Applying  the  covariance  sequence  model  of  (1.1)  to  the  co- 
variance  sequence  associated  with  the  square  root  of  a desired  spectral 
density,  one  may  desigr  recursive  digital  filters  with  a magnitude 
squared  frequency  response  that  approximates  the  desired  spectral 
density. 

In  Chapter  2 spectral  representations  are  given  for  stochastic 
processes  and  covariance  functions.  The  relations  between  the  proposed 
covariance  sequence  approximant  and  the  DFT  and  Pisarenko  decomposition 
are  explored.  A generalized  stochastic  process  synthesis  equation  is 
developed  for  a subset  of  those  processes  whose  covariance  sequences  are 
positive  real  linear  combinations  of  complex  exponentials. 

In  Chapter  3 we  develop  a modified  least  squares  fit  of  the  co- 
variance  sequence  approximant  of  (1.1)  to  a given  covariance  sequence. 

A decoupling  of  the  minimization  problem  for  the  complex  exponentials 
and  the  complex  weights  of  the  linear  combination  results.  An  efficient 
recursive  procedure  is  outlined  for  the  determination  of  a minimum  phase 
covariance  sequence  predictor  polynomial. 

In  Chapter  4 the  parametric  spectral  estimate  is  derived  that  is 
associated  with  the  covariance  sequence  approximant  of  this  thesis.  We 
show  that  harmonic  processes  in  white  noise  can  be  represented  as  a uni- 
form limit  of  this  spectral  estimate.  The  covariance  sequence  approxi- 
mant is  then  applied  to  finite,  known  covariance  sequences  as  well  as 
finite,  estimated  covariance  sequences  for  harmonic  processes  in  noise. 
As  a test  of  the  robustness  of  the  covariance  sequence  approximant  we 
apply  it  to  a Gaussian  covariance  sequence. 
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By  Parseval's  theorem  we  can  interpret  covariance  sequence  approxi- 
mation as  spectral  density  approximation,  so  that  covariance  sequence 
approximation  seems  the  natural  approach  to  designing  recursive  digital 
filters  from  magnitude  squared  frequency  responses.  An  efficient  design 
procedure  is  developed  in  Chapter  5.  Design  examples  are  given  for  some 
representative  lowpass  filter  designs  of  ARMA(2M,2M“)  type.  Important 
features  are  the  stability  of  the  filter  and  the  efficiency  of  the  de- 
sign algorithm. 

In  Chapter  6 conclusions  are  advanced  and  recommendations  for  ex- 
tensions are  proposed. 


t 
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2 SPECTRAL  REPRESENTATIONS  FOR  STOCHASTIC  PROCESSES 
AND  COVARIANCE  FUNCTIONS 

An  introductory  exposition  is  given  in  Section  2.1  for  general 
spectral  representations  of  wide  sense  stationary  stochastic  processes 
and  their  corresponding  covariance  functions  [K] , [Y],  [F],  [BR]. 

In  the  following  section  the  spectral  representation  for  the  co- 
variance  function  of  a linear  filter  driven  by  white  noise  is  treated. 
The  main  result  of  this  section  is  that  for  strictly  proper  ARMA(M,N) 
systems  the  spectral  representation  for  the  covariance  function  can 
always  be  written  as  a complex  linear  combination  of  complex  damped 
exponentials. 

Section  2.3  gives  the  spectral  representation  for  the  stochastic 
process  and  a stochastic  process  synthesis  equation  for  the  spectral 
analysis  procedure  based  on  the  Discrete  Fourier  Transform. 

The  spectral  representation  for  the  stochastic  process  and  the 
corresponding  stochastic  process  synthesis  equation  underlying  the 
Pisarenko  decomposition  spectral  analysis  procedure  are  treated  in 
Section  2.4. 

The  final  section  of  this  chapter  extends  the  synthesis  structure 
of  Sections  2.3  and  2.4  to  allow  the  generation  of  wide  sense  stationary 
processes  for  which  the  covariance  function  is  represented  by  a positive 
real  linear  combination  of  complex  damped  exponentials.  The  latter  is 
important  since  it  provides  a finite  parameter  representation  for  a 
large  class  of  wide  sense  stationary  stochastic  processes. 

2.1  Spectral  Representations  for  Weakly  Stationary  Processes 

A weakly  stationary  or  wide  sense  stationary  process  {X(t),tER} 


satisfies  the  following  moment  conditions: 


By  Parseval's  theorem  we  can  interpret  covariance  sequence  approxi 


approximation  seems  the  natural  approach  to  designing  recursive  digital 


filters  from  magnitude  squared  frequency  responses.  An  efficient  design 


procedure  is  developed  in  Chapter  5.  Design  examples  are  given  for  some 


representative  lowpass  filter  designs  of  ARMA(2M,2M“)  type.  Important 


features  are  the  stability  of  the  filter  and  the  efficiency  of  the  de 


In  Chapter  6 conclusions  are  advanced  and  recommendations  for  ex- 


tensions are  proposed 
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E X(t)  - E X(0)  ••  m ; -«  < t < °°  2.1a 

and 

E X(t+x)X(t)  « E X(t)X(0)  - Cx(t)  2.1b 

We  will  now  state  the  spectral  representation  theorem  for  wide  sense 
stationary  processes. 

Spectral  Representation  Theorem; 

Let  (X(t),  teR},  be  a real-valued  weakly  stationary  process  with 
E X(t)  » 0 and  covariance  function  Cx(x)  - E X(t+i)X(t).  We  further 
assume  that  Cx(x)  is  continuous  at  x-0,  in  which  case  Cx(x)  is  continu- 
ous for  all  x.  On  the  basis  of  these  assumptions  there  exists  a complex- 
valued random  set  function  Z(A) , called  the  random  spectral  measure  of 
the  process,  such  that  the  process  has  the  spectral  representation 
00 

X(t)  = / elxt  Z(dX)  2.2 

— oo 

Here  the  integral  is  a Riemann-Stleltjes  integral  and  the  sense  of  the 
equality  is  that  the  Riemann-Stleltjes  sum  converges  as  a limit  in  the 
mean  to  X(t) . 

The  random  spectral  measure  has  the  following  properties 

Z(-A)  - Z*(A)  2.3 

which  is  a consequence  of  assuming  X(t)  to  be  real-valued.  It  follows 
that 

E Z(A)  - 0 2.4 

since  E X(t)  -0.  If  F(A)  is  defined  by  F(A)  - E|z(A)|2,  then  F(A)  is 
a measure  called  spectral  distribution  and  the  following  relationship 
is  valid: 

E Z(A)Z*(B)  - F(A^B) 

Note  that  if  A and  B are  disjoint,  E Z(A)Z*(B)  - 0,  which  says 


2.5 
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that  the  random  variables  Z(A)  and  Z(B)  are  uncorrelated  (or  orthogonal). 
The  latter  two  expressions  can  be  put  in  the  following  operationally  con- 
venient form: 


E Z(dX)Z*(dy) 


Fx(dX)  if  y - X 


2.6 


l 0 if  y i X 

where  Fx(dX)  ■ F((-«»,X+dX])  - F((-«»,X]).  That  Is  dy  and  dX  are  thought 
of  as  being  Increments  about  y and  X such  that  if  y + X,  they  are  dis- 
joint. For  the  covariance  function  of  the  real-valued  wide  sense 
stationary,  zero  mean  process  we  then  derive  the  spectral  representation 
as  follows: 


Cx(t)  - E X(t+x)X(t)  - 

- E(  / eiX(t+T)Z(dX))(  / eiyt  Z(dy))* 

—oo  —oo 


= E 7 7 eiX(t+T)  e-iytZ(dX)Z*(dy) 

— OO  — OO 


■ / / e1^X'll>  teiXxE  Z(dX)Z*(dy ) 

— OO  —OO 

00 

= / eiXxFx(dX)  2.7 

—00 

This  result  is  known  as  the  Bochner-Khint chine  theorem. 

Discrete  time  parameter,  weakly  stationary  stochastic  processes 
have  a spectral  theory  which  exactly  parallels  that  given  for  continuous 
time  processes.  Let  (X(k),  k“0,±l,...}  denote  a weakly  stationary 
process  with  discrete  time  parameter.  Then  the  spectral  representation 
of  X(k)  is 

v.  ixk 

X(k)  - / e1AKZ(dX)  ; k«0,±l, . . . 


2.8 
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The  spectral  distribution  function  Is  given  In  terms  of  the  spectral 
distribution  F(A)  by 

FX(X)  - F((-tt,A])  2.9 

The  covariance  sequence  has  the  spectral  representation 
Cx(k)  - E X(j+k)X(j) 


iAk 


/ e““  Fx(dA) 


2.10 
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Note  that  the  only  real  difference  between  continuous  and  discrete 
time  parameter  stochastic  process  representations  is  the  Interval  of 


integration.  With  the  Lebesgue  decomposition  of  F (A)  we  can  give  a 

A 


more  intuitive  meaning  to  the  spectral  distribution  function  by  express- 


ing It  as  the  sum  of  a discrete  and  a continuous  component  [C] : 


F(A)  * Fd(A)  + FC(A) 


2.11 


The  continuous  component  can  be  further  decomposed 

FC(A)  - Fac (A)  + FSC(A)  2.12 

In  the  sequel  we  will  neglect  values  on  an  a>  set  of  probability  zero 
[D],  and  leave  the  singular  continuous  component  of  FC(A)  out  of  con- 


sideration. The  discrete  spectral  distribution  F°(A)  is  characterized 


by  the  spectral  function  {px(A^ ) > indicating  the  spectral  mass  concen- 


trated at  frequency  A^: 

Fd(A)  * E p (A.) 

j:  AjeA  A J 


2.13 


The  continuous  component  F (A)  is  determined  by  the  derivative  of  the 


spectral  distribution  function,  f_(A)  - FC(A),  where  fv(A)  is  called 


the  spectral  density  function.  Whenever  the  covariance  function  C (t) 

X 


is  absolutely  integrable  or  summable  there  is  corresponding  to  it  a 


continuous  spectral  density  function  f (A)  [D]. 

A 


I 
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In  case  the  spectrum  contains  both  discrete  and  continuous  compo- 
nents the  spectral  mass  or  power  in  a set  of  frequencies  A is: 

F(A)  - I PX(A  ) + / fx(A)dA  2.14 

j : A j eA  A 


For  this  mixed  spectrum  the  following  covariance  function  spectral 

representation  is  valid: 

IA.t 

Cx(t)  « Z e J Px(Xj)  + / e1AT  fx(A)dA 

Xj 

A similar  decomposition  of  the  random  spectral  measure  exists 


2.15 


Z(A)  - Zd(A)  + Zc(A)  2.16 

with 

E Z,(A)Z*(B)  =0  V A,B 
d c 

Defining  F(A)  as  before  this  leads  precisely  to  the  Lebesgue  decomposi- 
tion. Related  to  the  Lebesgue  decomposition  is  the  following  theorem 
due  to  Wold. 

Wold  Decomposition  Theorem: 

Let  {X(t),  t=0,±l, . . . },  be  a zero-mean,  weakly  stationary 
stochastic  process.  Then  X(t)  can  be  expressed  as  the  sum  of  two  zero- 
mean,  weakly  stationary  processes, 

X(t)  - U(t)  + V(t)  2.17 


such  that 

i)  the  process  U(t)  is  uncorrelated  with  the  process  V(t); 

ii)  U(t)  has  a one-sided  moving  average  representation 

00  OO 

U(t)  “la.  C(t-k)  with  a.-l  and  Z a?  < ® ; £(t)  is  white; 

k-0  K u 0 R 

iii)  the  V(t)  process  is  completely  determined  by  linear  functions 
of  its  past  values. 

The  process  V(t)  is  said  to  be  deterministic  and  the  process  U(t)  is 
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said  to  be  nondeterministic.  The  spectral  distribution  of  X(t)  can  be 
evaluated  in  terms  of  the  spectral  distributions  of  U(t)  and  V(t). 

With  (2.17)  and  property  i): 

Fx(dX)  - Fy(dX)  + Fy(dX)  2.18 

The  process  U(t)  has  continuous  spectral  distribution  due  to  property 
ii)  above  and  therefore  (2.18)  is  actually  the  Lebesgue  decomposition 
°f  FX^'  In  a11  cases  of  practical  interest  then,  the  deterministic 
component  is  simply  a stochastic  almost  periodic  function: 
n lX^t 


" XA  .U 

- Z Z e J j -«  < t < 
j— n 3 


where  *o»  *±1 ^±n  are  fixed  frequencies  and  the  are  complex 

random  variables.  For  V(t)  to  be  real-valued  and  weakly  stationary, 
we  need  the  following  properties: 

X.j  **  j=0, ±1, . . • ,±n  2.20a 

Z-j  ‘ ZJ  J-O.il.-.l"  2.20b 

I!  ZjZ*  - 0 for  j 1*  k 2 20c 

The  latter  results  in  the  following  spectral  representation  for  the 
covariance  function  of  a stochastic  almost  periodic  function 
n - iX  x 

C (t)  - Z E|Z  T e 3 2.21 

j=-n  3 

Using  the  properties  (2.20a)  aid  (2.20b)  this  can  be  written  as 
n 2 

C (t)  ■ Z 2o,  cosX.t  2.22 

j“0  3 3 

. (to) 

WhSre  Zj^  ’ rj  ((**)e  J in  (2.19).  This  means  that  V(t)  is  a linear 
combination  of  sinusoids  with  random  amplitudes  r^  (u>)  : (0,oZ)  and 
frequencies  Xj,  each  with  a random  phase  term  ^ (u)  : U[0,2n],  where 
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E r j (“>)  (u>)  - a2  6(j-Jl) 

E ^ (ai)^A(ai)  - 0 V 

E rjCwH^w)  - 0 V j,i  2.23 

By  virtue  of  the  Wold  decomposition  theorem  and  the  stochastic  almost 
periodic  function  of  (2.19)  we  can  now  consider  the  problem  of 
analyzing  sinusoids  in  noise  via  the  spectral  representation  of 
weakly  stationary  stochastic  processes. 

2.2  Spectral  Representation  for  Covariance  Functions  and  Linear 
Operations 

In  this  section  we  examine  the  effects  of  a linear  filter  on  the 

spectral  density  function  of  a so-called  white  noise  input.  A white 

noise  sequence  is  a sequence  of  uncorrelated,  zero-mean  random  variables 

2 

with  common  variance  a . This  sequence  is  a weakly  stationary  process 
with  covariance  function 


C(k)  = a26(k) 


6(k) 


1 k“0 

0 k#> 


2.24 


Since  obviously  (C(k)}ei^,  there  exists  a continuous  spectral  density 
function,  so  that  with  (2.10) 

C(k)  - / eiXk  f (A)dA  2.25 

-IT 

If  {C(k)}ei2*  as  In  most  applications,  then  the  numbers  C(k)/2ir  are 
the  Fourier  coefficients  of  the  Fourier  series  expansion  of  the 
function  f(A)  periodically  extended: 


f(A) 


2it 


£ 

k”-°° 


-iAk 

e 


C(k) 


2.26 
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Equations  (2.25)  and  (2.26)  are  the  discrete  time  version  of  the 
Wiener-Khintchlne  relations.  For  the  above  defined  white  noise  sequence, 
this  readily  yields 

f(A)  - o2/2tt  ; | X | <_  it  2.27 

Similarly  we  can  define  a generalized  weakly  stationary  process  in  the 

continuous  time  case  by  its  covariance  function 

C(t)  - o26(t)  2.28 

This  function  is  not  continuous  at  t_0,  but  has  a continuous  spectrum 

and  spectral  density 
2 

f (X)  ■ o /2u  ; -«»  < X < «°  2.29 

Denote  the  output  Y(t)  of  the  linear  filter  with  transfer  function 
H(X)  excited  with  input  X(t)  in  the  following  operator  notation 

Y(t)  = L{X(t) } 2,30 

The  action  of  a linear  filter  on  X(t)  is  completely  determined  by 
what  it  does  to  the  complex  exponential  function  eiXt  for  all  X,  which 
is  determined  by  the  transfer  function  H(X)  of  the  filter.  With  the 
spectral  representation  theorem  for  X(t) , we  can  therefore  write 
Y(t)  - L{X(t) } 

- / eiXt  H(X)  Zx(dX)  2.31 

The  latter  represents  the  output  'of  the  linear  filter  as  a weakly 
stationary  stochastic  process  with  spectral  measure 

Zy(dX)  - H(X)  Zx(dX)  2.32 

Thus 

Fy(dX)  * E |Zy(dX) | 2 

- |H(X)|2  Fx(dX)  2.33 

which  indicates  the  following  relations  between  the  spectral  functions 
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and  spectral  density  functions  at  the  Input  and  output  of  the  linear 
filter: 

PY(X)  - |H(X)|2  px(X)  2.34a 

fY(X)  - | H (X) | 2 fx(X)  2.34b 

Let  us  now  Investigate  the  covariance  sequence  structure  for  some 
well  known  processes,  resulting  from  white  noise  (0,o2)  driven  AR,  MA 
and  ARMA  systems  respectively. 

Assume  given  a stable  AR  system  with  transfer  function 

1 


1+a^z  ^+.  ..+8^  ** 


A(z) 


M -1 
n (l-p.z  l) 

1-1  1 


2.35a 

2.35b 

2.35c 


with  real  coefficients  a1  and  | p± | <1  (stability  requirement).  The 
spectral  density  function  for  the  output  then  satisfies  (2.34b)  and 
(2.27): 

2 


fAR<X>  “ 


1*0  2 


2n | A(e  A) | 


2.36a 


M 

2n  n (l-p.e“iX)(l-p  eU) 

i-1  1 1 


2.36b 


The  spectral  representation  for  the  output  covariance  sequence  (2.25) 

Is 


1 «1Xk  f. 


CAR(k)  " J fAB(X)dX  - 

»t  “ 


±. 

2»J  f M 


dz 


2.37 


n (*“?*) (1-p.z) 
1-1  1 1 


where  r is  a counter  clockwise  contour  lying  In  the  region  of  absolute 

convergence  of  the  Integral.  With  the  residue  theorem  we  then  find 

for  nonnegative  lags,  and  the  assumption  of  simple  poles: 

2 

C.p(k)  ■ a E residues  Inside  r 


For  negative  lags  use  the  definition  C.  (-k)  * CA  (k) . 

AK 

Assume  given  a moving  average  system  with  real  coefficients  b 
and  with  transfer  function 


When  the  above  MA  system  Is  driven  by  a white  noise  sequence  (0,o2) 
the  output  spectral  density  is  given  by 


Consequently  the  spectral  representation  for  the  corresponding  covariance 
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15 


N -1 

n (l-q  z X) 

- 1-1  1 

M -1 

n (l-p.z  L) 

l-i  1 


IpJ  < i 


2.A5c 


Driven  by  a white  noise  sequence  (0,q  ) the  resulting  output  spectral 
density  is  given  by 


f (x)  . qj  |MeU)|2 

ARMA(M,N)va'  2 it 


N 


2.46a 


2 H (l-q  e"1X)(l-q.eiX) 
_ q 1-1  3 J 


2ir  M -IX  IX 

n (l-p  e iA)(l-p.eiA) 

1-1  1 1 


2.46b 


The  spectral  representation  for  the  corresponding  covariance  sequence 
yields 

IT 


CARMA(M,N) (k)  “ /e  fARMk(M,N)(X)  dX 

“If 


2.47a 


_ _d_  i k-HM-H  jjvl 
2wJ  j,  M 


n (z-qjKl-qjZ) 


dz 


n (z-p.)(i-p.z) 
1-1  1 1 

For  k > 0 we  derive  via  the  residues  for  poles  inside  T: 

N 


2.47b 


M 


k-l+M-N  _ , 
z n (z-q  )(l-q  z) 

1-1  3 J 


- 0 * lln  (2'I>I)  IT- 

t-1  2^t  n (r-p.)  (l-p.z) 

1-1  1 1 


, 2 dn_1{znI(z)}  ,,  v 

+ a lim 1 * — I(n) 

z-K)  (n-1) ! dz 


2.48 


where  I(z)  in  the  Integrand  in  (2.47b)  and  n-N-hH-l-k  relates  to  a pole 
at  z-0  of  order  n.  The  indicator  function  I(n)  is  zero  except  when  n>0. 


Therefore 


m AT4*  (pr^Hi-qjP,) 

CARMA(M,N) (k)  “ f M p* 

(1-p?  (prpi)(1_pip£) 


+ lim 


.v-1  A (z'qj)(1-v>| 


z-»0  (n-l)!dz 


n-1  1 M 


I(n)  2.49 


n (z-p  )(l-p  z) 
i*=l  1 1 


The  covariance  sequence  symmetry  property  again  defines  the  covariance 
sequence  for  negative  lags.  From  the  above  results  we  deduce  that  the 
covariance  sequence  for  an  ARMA(M,N)  process  with  simple  poles  can  be 
represented  in  the  form  of  a complex  linear  combination  of  complex 
exponentials  whenever  the  numerator  order  N is  at  most  one  less  than 
the  denominator  order  M,  which  we  will  denote  M in  the  sequel.  For 
these  so-called  strictly  proper  rational  transfer  functions  the  follow- 
ing covariance  sequence  model  is  valid: 


CARMA(k)  " 1 \ PZ  5 V k 2.50 

i-1 

M N 

Because  the  coefficients  {a^}^  and  {5^}^  are  real,  we  have  the  property 
M 

that  are  either  real  or  occur  in  complex  conjugate  pairs. 

At  this  point  we  note  that  in  particular  the  covariance  sequence 
expressions  (2.39)  and  (2.49)  were  derived  to  show  when  expression 
(2.50)  would  be  valid.  If  one  desires  to  compute  numerically  the 
covariance  sequence  for  a general  ARMA(M,N)  system  driven  by  white 
noise,  it  is  much  more  convenient  to  solve  the  linear  set  of  general- 
ized normal  equations  of  order  M,  extend  the  resulting  sequence  to 
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order  max(M,N)  and  subsequently  use  the  autoregressive  recursion  of 
the  denominator  to  compute  as  many  of  the  positive  going  lags  as  desired. 
This  approach  Is  discussed  In  [DBE]. 

2*3  The  DFT  and  its  Synthesis  Process  Spectral  Representation 

Assume  we  have  the  following  partial  sequence  realization  x(0), 
x(l) , . . . ,x(N-l) . For  the  periodic  extension  of  this  partial  realiza- 
tion write  the  Fourier  series  expansion  as  follows 
N-l  IX  k 

X(X  ) - Z x(k)e  n 2.51a 

n k=0 


where 


. N-l  -IX  k 

x(k)  = - Z X(X  ) e n 2.51b 

n=0  n 


with 


X - 


2im 


0 < n < N-l 


n N 

These  relations  constitute  the  Discrete  Fourier  Transform  (DFT)  pair, 

which  can  be  implemented  very  efficiently  on  a digital  machine,  due 

iXQk  iXN_nk 

to  the  symmetry  properties  of  e and  e “ . The  algorithms 

referred  to  are  known  as  FFT  algorithms.  Recall  that  the  spectral 
representation  for  a stochastic  almost  periodic  function  was  given 
in  (2.19): 

IX. t 


V(t) 


n 

Z 


j«-n 


7 e ^ 

j 


— °o  < t < 00 


Due  to  the  symmetry  property  of  the  {X ^ } and  the  orthogonality  property 
of  the  {Zj},  as  in  (2.20),  we  recognize  that  x(k)  in  (2.51b)  is  a 
particular  realization  of  the  stochastic  almost  periodic  process  with 

harmonically  related  frequencies  X..  The  power  at  frequency  X is 

2 J ^ 
reflected  in  E{ | X(X ^ ) /N | }.  For  a particular  partial  realization  of 
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a process,  the  DFT  therefore  merely  yields  an  estimate  of  the  power  at 
the  harmonically  related  frequencies  X ^ . 

It  is  worthwhile  to  note  that  the  covariance  sequence  corresponding 
to  the  analysis  model  underlying  the  DFT  is  represented  by  a positive 
real  linear  combination  of  harmonically  related  complex  exponentials  as 
in  (2.22).  DFT  spectral  analysis  or  the  perlodogram  is  therefore 
recognized  as  a special  case  of  the  model  presented  in  (2.50). 

2.4  Pisarenko  Decomposition  and  its  Synthesis  Process  Spectral 

Representation 

Let  us  assume  a very  practical  situation  where  we  have  a stochastic 
almost  periodic  process  X(k)  in  additive  noise  N(k),  that  is  orthogonal 
to  the  process  X(k) . The  covariance  sequence  is  now  given  by 
c(k)  - + Vk)  (X*  + N*)) 

' E{lw0  + E(w‘i  + + E{N(+kV 

■ cx(k)  + CB(k)  2.52 

C^tk)  is  given  by  (2.22)  and  C (k)  depends  on  the  dynamic  properties 
of  the  noise  generating  process.  For  the  particular  case  where  the 
assumption  is  made  that  N(k)  represents  a white  noise  process,  the 
above  results  in  the  model  underlying  the  Pisarenko  decomposition  [P]. 

2 n 

C (k)  « o„  6(k)  + E p,  cosX  k 2.53 

P w j-0  J 3 

The  Pisarenko  decomposition  is  thus  seen  to  decompose  a finite  covari- 
ance sequence  realization  into  a white  noise  component  and  a positive 
real  linear  combination  of  complex  exponentials  which  are  not  necessarily 
harmonically  related.  If  the  underlying  covariance  sequence  generator, 
or  the  process,  satisfies  the  assumed  model,  then  the  parameters  in 
(2.53)  can  be  determined  uniquely  according  to  a theorem  by  Carathdodory 
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[GS] . The  stochastic  process  synthesis  equation  can  obviously  be 
written  by  Inspection  as  follows: 

n 1A  ,k 


X (k,o>)  - E Z (w)e  i + N(k,w) 
P j--n  J 


2.54 


with 


I Pr 


EZ, 


Pj/2 


j=0 

j*> 


E N 


'N 


Referring  again  to  the  model  In  (2.50),  we  note  that  a real  pole  p 

£ 

can  approximate  arbitrarily  closely  the  white  noise  Imposed  6(k)  when 
the  pole  radius  goes  to  zero.  The  model  underlying  the  Pisarenko 
decomposition  then  becomes  a special  case  of  the  model  presented  in 
(2.50)  as  well. 

2.5  A Generalized  Synthesis  Equation  for  Wide  Sense  Stationary  Processes 
Now  the  Interesting  question  arises  as  to  whether  the  stochastic 
process  synthesis  structures  as  Indicated  for  periodogram  analysis  and 
Pisarenko  decomposition,  can  be  extended  so  as  to  allow  damped  complex 
exponentials  in  the  spectral  representation.  Let  us  therefore  investi- 
gate the  following  stochastic  process  synthesis  equations: 


M j*>)  jX  k 

X(k,u)  - E v . (u>)e  m e “ 
m=l  m’k 


2uf  T 
m 


2.55a 


Here  the  complex  random  variable  v . (w)  is  orthogonal  to  the  random 

m,K 

phase  variable  $m(u>) , and  satisfies  a first  order  autoregression  [S]: 

"..t'"1  ■ + £»,k  2-55b 


We  furthermore  choose 


1-1 6, 


m1 


6(m-t)6(k-j) 


2.55c 
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Let  us  determine  the  covariance  properties  for  v . (u),  under  the 

m9ic 

assumption  that  v n(w)  is  orthogonal  to  e.  , for  all  £,k: 

E{Vk+j(“)Vm,k(u)}  ’ 

k+1  k+j-1  * ^-1  *o  * 

- E{(0kt;,v  n(u>)  + Z py  e ...  )(g  v rt(u>)  + Z 8 p e , „)} 

m n-0  m m m»^+j-n  m m,0  ' 9mFi m Mm  m.k-i'' 


£-0 


o(oi)vV  0(u)}  + E * bV^Ip  |2E  {e  . e*  #} 
m m m,0  m,0'  . n n m m 1 m1  m,k+j-n  m,k-£ 


£-0  n-0 
k-1  k+j-1 


8>jX,o(*>V»«)  + Aedl6ml2t|>'„l(1-I6j2) 


1-18.1 


2k 


ei|sJ2kE{v  (m)v*  (u>)}  + eilum|(l-|Bm|2)  

m m'  m,U  m,0  m1  m'  1 m'  . i i2 

' m 

8>J  + fil»«l2k[Wv.,oWvIt0M>  - l».U 


2.56 


Covariance  sequence  stationarity  for  v , (u)  is  seen  to  be  obtained 

m,K 

with 


E{vm,0(w)vm,0(“)}  " |wml 

From  (2.56)  and  wide  sense  stationarity  of  v . (u)  we  note: 

m,k 

E{vm,k+j(“)vt,k(t0)}  " 

The  autoregression  of  (2.55b)  furthermore  yields: 


k-1 


- E{gk\>  (u>)v  (o>)  + Z 8^  p e . . (w)  } 

m m,£  m,£  m m m,k-j+£  m,£v  ' 

■ 3*  E{v  „(w)v*  (w)} 

m m,  £ m,  £ 

■ < K\ 


2.57 


2.58 


2.59 
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For  the  covariance  sequence  of  process  generator  (2.55a)  we  finally 
derive : 


C(k)  - E{XU+k,w)X*U,u>)} 

M J ♦_(«))  jXma-Hc)  M * -j*,(u)  -jX.i 

* E{  E v j.  (u>)e  m e m Z v.  tf(m)e  j e i } 

m«l  ’*■  j-1 


M 

■ Z E{v 
m“l 


* JX  k 


M k JAn.k 

1 m1  m 

m-1 


M 

* |U 

m*l 


m"eml  e 


k J(Xm+ar*em)k 


2.60 


j $„(<*>)  jXmk 

e ) real  or  In  conplex  conjugate  pairs  the 


With  {v  . (uj)e  m 
m,K 

process  X(k,u»)  is  real.  Note  from  (2.60)  that  the  linear  weight 
involved  is  positive  real,  and  with  the  generator  equations  (2.55)  we 
can  therefore  synthesize  a weakly  stationary  stochastic  process  in 
the  class  of  ARMA.(M,N)  processes,  which  exhibits  covariance  sequences 
that  are  positive  real  linear  combinations  of  damped  complex  exponen- 
tials. This  point  may  be  clarified  as  follows.  Consider  an  ARMA(M,M~) 
system  with  real  coefficients  and  simple  poles,  so  that  its  impulse 
response  may  be  written: 


M 


"k  ■ A Bj 


V k > 0 


2.61 


For  the  covariance  function  of  the  stable  system  we  derive: 


ck  - ,f0  h*  Vlkl 


M M “a  A+lkl 
z b i b £ pJ'p:  |k| 

i-1  1 j-1  J A-0  1 J 

**  **  1 Ijtl 

£ [B,  E B,  T ■ ■ ■ ] P,|K| 

j-1  j i=l  1 1-*iPj  J 


M lBll  M Bi  Ikl 

£ [ o + B,  £ ■---■*  ] P,11 

J-1  l-|p  |2  J 1-1  1_PiPJ  J 

J 


2.62a 


2.62b 


where  t*  denotes  the  exclusion  of  1 such  that  B^B^ . Note  that  by 
defining  the  term  In  square  brackets  of  (2.62a)  to  equal  the  general 
covariance  sequence  model  of  (2.50)  for  strictly  proper  ARMA(M,M-) 
systems  has  been  derived  In  an  alternative  way.  The  Important  inter- 
pretation of  (2.62b)  Is  that  the  second  term  In  square  brackets 
represents  Interactions  between  the  different  modes  of  the  ARMA(M,M~) 
system,  resulting  In  a linear  weight  coefficient  that  is  not  necessarily 
positive  real.  Thus  these  interactions  may  not  be  represented  in  the 

j* 

synthesis  model  of  (2.60).  For  a real  process  (Bo  } are  either  real 

m 

or  occur  in  complex  conjugate  pairs,  and  (2.60)  thereby  represents  a 
positive  real  linear  combination  of  nonnegative  definite  covariance 
functions  [Y],  which  implies  that  C(k)  of  (2.60)  is  nonnegative  definite 
[W].  By  the  Herglotz  theorem  it  is  only  for  nonnegative  definite 
functions  that  the  spectral  representation  of  (2.10)  exists  [ BU ] . 

From  (2.62b)  we  see  that  (2.60)  does  not  represent  the  most  general 
linear  combination  of  damped  exponentials  that  represents  a nonnegative 
definite  covariance  function.  A further  generalization  would  be  achieved 
when  |pjn|in  (2.60)  were  replaced  by  real  or  complex  conjugate  pairs  of 
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linear  weights  such  as  indicated  in  (2.50)  derived  from  strictly 
proper  ARMA(M,N)  systems  driven  by  white  noise.  The  gain  derived  from 
the  more  general  representation  of  (2.50)  is  illustrated  by  the  spectral 
representation  theorem  for  wide  sense  stationary  processes  as  given  by 
Yaglom  [Y]: 

Given  any  wide  sense  stationary  process  £(t),  any  e > 0 and 

any  T > 0,  there  exist  random  variables  C^, 

are  uncorrelated  in  pairs  (the  number  n depends  on  c and  T, 

of  course)  and  real  nunfters  A,, A A such  that 

x z n 

n iA,  t „ 

E|  5(t)  - £ L e \ < e 2.63 

k-1  x 

for  any  t e[-T,T]. 


1 * * • i ^ which 


Every  wide  sense  stationary  stochastic  proqess  can  be  approximated 
arbitrarily  closely  in  mean  square  sense,  by  a stochastic  almost 


periodic  function,  but  to  achieve  better  accuracy,  in  general  the 
number  of  harmonic  components  has  to  be  increased  and  the  differences 
A^+^  - A^  between  neighboring  frequencies  has  to  be  decreased. 

For  covariance  sequences  as  derived  above  for  strictly  proper 
ARMA(M,N)  systems  Increased  accuracy  would  indeed  force  n-*»  and  A^.  - 
A^  -*■  0 in  (2.63).  This  would  leave  us  with  an  infinite  parameter 
representation  for  a clearly  finite  parameter  problem  indicated  by  (2.50). 
It  is  exactly  in  this  finite  parameter  spectral  representation  that  we 


find  a powerful  tool  for  parametric  spectrum  analysis  (Chapter  4)  and 
recursive  digital  filter  design  (Chapter  5). 

2.6  Diagram  for  Generalized  Stochastic  Process  Synthesis 

Generation  of  the  process  X(k,u>)  of  (2.55)  may  be  achieved  by  the 
implementation  in  Figure  2.1.  The  random  variables  can  be  considered 
real  and  are  orthogonal: 


V(u)  1L  em,k(u)  1L  vm,0(u) 

Their  respective  statistics  are  given  by: 

♦ (w)  - U[0,2ir] 
in 

i-lej2 

em,k(u)  * (0*  ■_,|yniT) 

v»,o(u)  ~ <°»kl> 


2.65a 


2.65b 


2.65c 


The  independent  branches  each  generate  a real  process,  or  they 
occur  in  independent  complex  conjugate  pairs  so  that  the  pair  generates 
a real  process.  The  sum  of  all  the  branches  is  then  real.  The  param- 
eter choices  for  the  particular  process  synthesis  equations  are  as 
follows: 

For  DFT  related  process  synthesis: 

an  6-1 

m 

, 2 inn  _ „ , 

Xmm  T~  ; 0 - ® - N_1  2.66 

For  Pisarenko  decomposition  related  process  synthesis: 

B - 1 V m + 0 

in 

Xm  not  nessarlly  harmonically  related 
white  noise  branch  for  m ■ 0: 

x0-o 

♦0  (o>)  - 0 

Bq  - 0 (pole  at  z«0) 
v0,0(u)  = 0 

W0  real  2.67 

The  generalized  process  synthesis  is  indicated  in  Figure  2.1.  Note 

that  if  (X  } are  to  indicate  the  constituent  frequencies,  (B  ) should 
“i  in 

be  taken  real  in  (2.60),  which  obviously  does  not  form  a restriction. 


■ 
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3 A MODIFIED  LEAST  SQUARES  PROBLEM  AND  ITS  APPLICATION  TO 
COVARIANCE  SEQUENCE  APPROXIMATION 
First  we  pose  the  general  least  squares  problem  for  approximating 
covariance  sequences  and  discuss  several  ways  to  solve  or  avoid  it. 

In  Section  3.2  the  approach  taken  in  this  thesis  is  outlined. 

First  the  complex  exponentials  for  the  proposed  covatiance  sequence 
model  are  found  by  minimizing  the  covariance  sequence  linear  prediction 
errors.  Subsequently  the  corresponding  linear  combination  weights  are 
determined  in  a general  least  squares  fashion.  This  necessarily  sub- 
optimal  approach  leads  to  linear  systems  of  equations  and  a polynomial 
rootfinding  problem  to  be  solved.  The  somewhat  ad  hoc  nature  of  this 


approach  precludes  the  derivation  of  a simple  spectral  error  criterion. 

A recursive  solution  to  the  covariance  sequence  predictor  is 
derived  in  Section  3.3  and  its  desirable  properties  are  enumerated  in 
the  final  section. 

3.1  General  and  Modified  Least  Squares  Problems 

For  our  covariance  sequence  model  we  choose  the  complex  linear 
combination  of  damped  complex  exponentials  found  in  the  previous 
chapter: 

M Ikl 

R.  - I A.p . Vk  3.1 

i-1  1 1 

Clearly,  R_^  - R^,  and  as  derived  before  {A^.p^}  are  either  real  or 
occur  in  complex  conjugate  pairs.  There  are  several  reasons  one  might  *\ 

\ 

choose  this  model: 

(1)  Such  a covariance  sequence  model  arises  naturally  from  a 
strictly  proper  ARMA(M,N)  process  driven  by  a white  sequence. 

. A 
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(2)  This  model  preserves  second  order  statistical  properties  in 
converting  from  continuous  time  to  discrete  time.  A 
strictly  proper  continuous  time  process  described  by  an  Mtk 
order  differential  equation,  transforms  under  those  preserva- 
tion conditions  to  a discrete  time  system  with  M poles  and  at 
most  M-l  zeros  [PS], 

(3)  For  physically  realizable  systems  it  seems  natural  to  repre- 
sent the  covariance  sequence  by  a complex  linear  combination 
of  complex  exponentials. 

(4)  The  model  corresponds  to  a generalized  stochastic  process 
synthesis  structure  based  on  positive  real  linear  combinations 
of  damped  complex  exponentials. 

In  fitting  the  model  of  (3.1)  in  a least  squares  sense  to  a given 
symmetric  sequence  (C^)*  k— 0,  ±1,  ±2,...,±L,  a simultaneous  minimiza- 
tion with  respect  to  and  p^  has  to  be  performed  on  the  following 
error  criterion: 

L 2 

e - z {C.-R.  r 

k— L K K 

L M Ikl  2 

“ I (C.  - E } 3.2 

It— L i-1  1 1 

Equating  the  partial  derivatives  of  E to  zero  yields: 


0mitm  * {V*  AipJk|}pi*k* 

3Aj  k— L k i-1  1 1 J 


; j-l,2,...,M  3.3a 


■jj~  “ E <V  l A p Jkl}A  |k|pjkl"1  ; j-1,2 M 

Pj  k=-L  i-1  1 1 J J 


The  above  equations  are  the  discrete  time  versions  of  the  Aigrain- 
Willlams  equations  [AW],  with  the  distinction  that  they  are  derived 
here  for  approximating  the  two-sided  covariance  sequence  {C^}  instead 
of  a one-sided  causal  impulse  response  sequence  {h^} . Hote  that  equa- 
tions  (3.3.)  .re  linear  la  (*t>“  For  my  give,  .et  of  pole.  (p^J 
the  corresponding  llne.r  combination  might.  (AjlJ  c.n  be  determined  lu 
straightforward  linear  fashion  by  solving  (3.3a): 
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and  write  (3.3a)  as: 

T T 

P Q C = P Q P A 


3.5b 


3.6 


Then  the  solution  for  the  linear  weights  is 

A = [PTQP]_1PT  Q C 3.7 

which  is  recognized  as  a weighted  least  squares  fit  for  covariance 
lags  k ^ 0.  The  weighting  is  entirely  due  to  the  two-sided  approxima- 
tion of  a real  symmetric  sequence  with  an  approximating,  real  svnmetric 
sequence.  The  approximation  for  turns  up  once,  whereas  the  approxi- 

mation for  C^,  k > 0,  turns  up  once  again  in  the  approximation  of  C 
resulting  obviously  in  the  weighting  matrix  Q above.  The  matrix  P is 
a Vandermonde  matrix  and  with  the  assumption  that  all  roots  {p^}^  are 

distinct  P is  full  rank  [C].  With  the  given  definition  for  Q the 
T 

determinant  of  P QP  differs  from  zero  and  a unique  solution  for  A 
results  from  (3.7)  [ND]. 

The  equations  (3.3b)  can  be  written  in  a similarly  reduced  form: 


'PTQC-'PTQPA  3.8 

Here  the  definitions  for  Q,  C,  P and  A remain  as  before,  and  'P  is 
defined  as 


M 

E 

j'l 


3.9 


Thus 
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M 3PT 
Z A.  Q C 

j-1  J 3PJ  " 


M apT 

Z A Q P A 
-1  J 


3.10 


J-1 


This  expression  is  nonlinear  in  {p^}.  In  the  very  special  case  where 
the  sequence  {Cfc}  has  an  exact  representation  of  the  form  (3.1)  the 


solution  obviously  satisfies: 


C = P A 3.11 

In  the  following  section  we  will  show  how  this  solution  can  be  obtained 
as  a special  case  when  the  model  exactly  matches  the  desired  covariance 
sequence  {C^}  over  the  entire  Interval  of  observation.  The  latter  is 
the  case  when  satisfies  an  ABMA(M,M-)  model,  or  when  the  chosen 
model  order  M is  so  large  that  an  exact  match  of  is  obtained  over 
the  interval  of  approximation. 

The  nonlinear  minimization  problem  outlined  above  generally 
requires  iterative  solution  methods,  several  of  which  have  been 
suggested.  Based  on  geometric  orthogonality  concepts,  an  orthonormal 
basis  can  be  generated  Iteratively,  and  the  projection  theorem  would 
lead  to  the  best  approximation  [MH] . Another  iterative  procedure  arises 
from  a system  identification  point  of  view  [SMC],  where  an  estimate 
(Aj_i(z))  for  the  denominator  polynomial  filter  is  used  to  prefilter 
in-  and  output  of  the  system  to  be  identified.  The  next  estimates  for 
denominator  (A^(z))  and  numerator  (Q^(z))  polynomials  are  obtained 
simultaneously  using  either  of  the  following  equivalent  generalized 
least  squares  formulations: 


min  ^ 
Ai’Qi 


Q±(z) 


A^z) 


'Ai  i(z)  ” A^Cz) 


H(z) 


2 dz 
z 


min 

Ai,Q 


i 


Qa(z) 

Ui(z) 


- H(z) 


2 1 A1(Z>  1 2 dz 
Ai-l(z) 


3.12a 


z 


3.12b 


\ 
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The  new  denominator  polynomial  filter  Is  used  to  Iterate  the  procedure. 
Experimental  results  have  shown  that  this  process  often  converges  quick- 
ly but  also  that  the  danger  of  divergence  is  imminent  when  too  many 
parameters  are  specified  [ST].  The  latter  may  explain  why  a general 
convergence  proof  has  not  been  given  since  appearance  of  the  original 
paper.  A similar  linear  iterative  process  with  improved  convergence 
properties  has  been  suggested  [Ml] , and  requires  only  one  mode  of 
iteration.  The  relations  between  [MH]  and  [SMC]  are  described  in  [M2]. 

The  criterion  in  (3.12)  is  based  on  a suggestion  by  Kalman  [K] 
and  is  known  loosely  as  the  modified  least  squares  criterion.  A recent 
approach  capitalizes  on  the  resulting  tractable  minimization  problem  to 
find  an  approximating  ARMA(M,N)  system  based  on  finite  impulse  response 
and  covariance  sequences  [MR].  In  the  modified  problem  one  seeks  to 
minimize  the  following  error  with  respect  to  denominator  coefficients 
{a^}  and  numerator  coefficients  {q^  > : 

E =*  / |H(ej“)A(ej“)-Q(eJui)|2  dm  3.13 

-Tf 

Note  that  the  important  difference  with  the  familiar  least  squares 
formulation  is  the  introduction  of  a spectral  weighting  function 
|A(e^U)|2,  which  happens  to  be  the  magnitude  squared  response  of  the 
yet  unknown  denominator  polynomial.  The  resulting  error  criterion 
represents  a quadratic  minimization  problem.  It  leads  to  a condition 
for  the  numerator  coefficients  depending  on  impulse  response  data  and 
the  denominator  coefficients.  When  this  condition  on  the  numerator 
coefficients  is  substituted  in  the  error  expression,  it  results  in  a 
constrained  quadratic  minimization  problem  for  the  denominator  coeffi- 
cients only.  The  authors  give  an  efficient  algorithm  of  Levinson  type 
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to  find  the  filter  coefficients.  The  most  Important  requirement  for 
the  above  procedure  to  yield  a stable  solution.  Is  that  the  finite 
sequences  {h^)  and  {C^}  be  consistent.  One  way  to  take  care  of  this 
requirement  [SL]  is  to  find  a very  high  order  autoregressive  approxima- 
tion first,  using  say  the  autocorrelation  method  of  linear  prediction 
[MA] . The  resulting  high  order  system  is  used  to  generate  consistent 
impulse  and  covariance  sequences,  before  applying  the  Mullis-Roberts 
algorithm  to  yield  stable,  low  order  ARMA(M,N)  systems. 

Time  domain  approaches  to  pole-zero  modeling  have  been  around  for 
a long  time.  Prony  devised  a method  for  functional  interpolation  [PR] , 
which  was  extended  later  to  include  least  squares  fits  [H],  The  original 
method  by  Prony  was  shown  to  be  equivalent  to  Padd  approximation  [WM] . 
Later  time  domain  approaches  have  been  termed  rediscoveries  of  Prony’ s 
method  [SH],  [BP].  These ' approaches  are  all  based  on  the  fact  that  the 
tail  of  the  impulse  response  of  an  ARMA(M,N)  system  satisfies  a pure 
autoregression.  The  latter  property  allows  one  to  find  the  system  poles 
first,  and  thereafter  the  zeros  for  given  poles.  These  procedures  have 
been  very  popular  in  spite  of  the  fact  that  in  general  stability  cannot 
be  guaranteed  for  the  resulting  filters.  The  procedure  in  the  following 
section  is  quite  similar  to  the  above  time  domain  least  squares 
approaches.  Differences  lie  in  the  application  to  two-sided  covariance 
sequence  approximation  and  the  possibility  to  carefully  monitor  the 
stability  of  the  recursively  developing  denominator  polynomial,  so  that 
stability  can  be  guaranteed. 

3.2  A Covariance  Sequence  Prediction  Approach 

The  question  arises  as  to  whether  the  general  nonlinear  least 
squares  problem  can  be  approximated  as  far  as  the  error  criterion  goes. 
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but  avoided  as  far  as  iterative  solution  methods  are  concerned.  Remember 
that  the  covariance  sequence  model  arose  naturally  from  an  AKMA(M,M~) 
system  driven  by  a white  sequence.  We  can  therefore  write  the  differ- 
ence equation  for  such  a system: 


1 ai  yk-i  " z bi  Vi  ; a0  " 1 

i“0  1 j=0  J k J 0 

The  impulse  response  then  satisfies  the  same  difference  equation 


1 ai  Vi  “ E bi  Vi 

i=0  1 1 j-0  i k 3 

and  with  the  definition 


Ck  = 1 hihi+k  3.16 

i=0  1 1 K 

the  covariance  sequence  for  an  AKMA(M,M”)  filter  is  found  [KPSS]  to 
satisfy: 

M M-l-k 

lVi°k-i'  slo  "iV*  ' ki°  3-17 

It  is  Important  to  note  that  the  covariance  sequence  of  this  ARMA(M,M“) 
system  satisfies  a purely  autoregressive  difference  equation  for  lags 
k >_  M.  We  will  use  this  property  in  order  to  find  the  denominator 
polynomial  coefficients  {a^}^,  *rom  which  the  poles  {p^}^  can  be  found 
by  polynomial  rootfinding.  From  (3.5)  note  that  if  the  complex  expo- 
nentials are  distinct,  the  (A.}?  can  be  found  such  that  Cn,..  C 

xx  0 ’ j^-l  dre 

exactly  matched.  Under  this  condition  the  minimization  problem  of  (3.2) 
gets  replaced  by 


min  I {c,  — R 
{p^  k=M  k * 
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A weighting  of  these  approximation  errors  by  the  autoregression  of  the 
polynomial  corresponding  to  {p^}  analogous  to  (3.13)  leads  to 
L M M 

min  E { E a C.  . - E a.  R.  Y 3.19 

{a^  k-M  i-0  1 i-0  1 

Recognize  that  the  approximating  sequence  {R^J  exactly  satisfies  the 
autoregressive  recursion  as  in  (3.17).  An  alternative  interpretation 
of  this  procedure  is  therefore  that  we  are  looking  for  the  best  linear 
predictor  for  the  covariance  sequence: 


min  E (C.  + E a,  C,  . } • 

{aj  k-M  k 1-1  1 *-*• 

Taking  derivatives  with  respect  to  (a.)  leads  to  the  following  systc 

3* 

of  linear  equations: 


{Ck  + ai  Ck-i}  Ck-j"°  ; J"1*2* •••.** 


Let  us  write  them  out: 


St-1  CM  CL-1  Sl-1  CM-2*  C0 

Si- 2 . Si  Sl-1  C1 


c0  c1  --.c 


S-l  S.-2*  ’ • S^M 


Define: 


Sl-1  Sl-2*  C0 
Si 


St-1***  CL-1 


c0  ... 


3.22a 


S-l  CL-2*  CL-H] 
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a - [a,,  a2,  ...  .a^ 

-T  “ [CH’  Sn-1  • • • • tCj  3.22b 

which  yields  the  following  matrix  equation 
T T 

C C a = -C  c 3.23 

T 

The  elements  of  C C are  covariance  estimates  of  the  covariance  sequence, 
given  by 


{(TC).  - 1 C,  C.  . ; l<i,J<M 

13  k-M  3 


and  it  is  easily  proved  that  CTC  is  nonnegative  definite.  Of  course 
(3.23)  can  be  solved  by  general  linear  systems  routines,  but  its 
symmetry  can  be  used  to  advantage  and  the  more  efficient  Choleski  de- 
composition performed.  The  linear  system  of  equations  (3.23)  has  a 

T 

unique  solution  when  the  determinant  of  C C does  not  equal  zerc,  a 

condition  met  if  and  only  if  C is  full  rank.  In  case  the  underlying 

system  for  has  order  less  than  M,  the  leftmost  column  of  C can  be 

expressed  as  a linear  combination  of  the  other  columns  and  C is  not 

full  rank.  Many  solutions  would  be  possible  in  that  case.  We  note 

here  that  in  practical  cases  {C^}  does  not  exactly  satisfy  an  Mth  order 

difference  equation  and  (3.23)  will  have  a unique  least  squares  solution. 

T 

Furthermore,  even  when  C C is  singular,  due  to  an  over  estimation  of  the 

underlying  autoregressive  system  order,  the  event  can  be  effectively 

anticipated  by  a scheme  that  recursively  updates  the  autoregressive 

coefficients  and  the  system  order.  By  keeping  track  of  the  covariance 

T 

sequence  prediction  error  at  each  step,  singularity  of  C C can  be 
detected.  We  point  out  here  that  (3.23)  represents  the  covariance 
method  of  linear  prediction  for  predicting  the  tall  (k>M)  of  an  other- 
wise synmetric  real  covariance  sequence. 
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The  covariance  method  of  linear  prediction  can  be  derived  under 
the  assumption  that  the  portion  of  the  signal  from  which  the  predictor 
parameters  are  computed  Is  nonstationary  [MM].  Let  (Rc(£,k)}  denote  the 
nonstatlonary  covariance  function  for  the  covariance  sequence  {C^}^. 

Let  {R  (i,k)  } denote  the  covariance  function  for  the  error  sequence 
{efc}  observed  when  the  linear  prediction  filter  A(z)  is  excited  by  {Cfc}. 
Thus  Rc(i,k)  may  be  written 


R U,k) / / Q («,«‘')e^  <ft,*""w>lt)d(i»d<»»' 

C (2tt)Z  -v  -ir  c 


3.25 


with 


Q (w,a)')  - Z Z R (l,k) 
k—  A— 

It  follows  easily  [PA]  that 

Q. <«.»')  - Qc(m,u.')A(eJ<i,)A*(eJ“') 

Thus  the  error  variance  E » Re(0,0)  may  be  written 


w * 


E - 


■J  / / Qc(w,cu')A(e-*“)A*(e^U)  )du>dw' 


(2w)  -v  -w 

With  substitution  of  the  linear  prediction  filter 


3.26 


3.27 


3.28 


B “ — ^~2  / / Qc(w,uO[l  + Z a.  e“j“k][l  + I a.ejw^ 
(2v)z  -n  -»  c k-1  ^ £-1  1 


]do)dw' 


3.29 


Minimizing  E with  respect  to  the  predictor  coefficients  yields: 

IT  n M ^ 

/ / Qc(“»«')  lw{L+  Z a eiU>  H*Ji“  (1+  Z a.  e“jwk}]do)du)'-0; 
-*  “»  t-1  k-1  ^ „ 


1-1, 2,..., M 
« 3.30 


The  definition  of  (3.26)  allows  us  to  write 
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M 

£ \ R (-i,-k)  = -R  (-1,0)  ; 1 < i < M 3.31 

k«l  * c c “ “ 

Equations  (3.31)  are  known  as  the  generalized  normal  equations. 

A 

With  Qc(u),oO  chosen  to  be 


Qc(u,o>') 


A(eJ“)A*(eJ“') 


3.32 


the  spectral  error  criterion  can  be  rewritten  as 

e - & 7 J 


dwdo)' 


-ir  -ir  Qc(o>,<iO 


3.33 


Minimization  with  respect  to  the  predictor  coefficients  of  either 
(3.29)  or  (3.33)  is  seen  to  yield  the  solution  of  (3.31).  Similarities 
are  present  here  with  the  analysls-by-synthesis  procedures  when  approxi- 
mating with  a 2D  all-pole  model.  Note  that  if  Qc(*»*)  is  larger  than 


Qc(’»*)  a significant  error  contribution  results.  If  on  the  other  hand, 
Qc(*»  *)  Is  smaller  than  Qc(*,*)  a relatively  small  enror  contribution 
takes  place.  As  a result,  the  best  approximation  will  occur  where 
Qc(*» •)  is  high.  Another  observation  to  be  made  relates  to  the  unifor- 
mity of  the  error  criterion.  Due  to  the  ratio  of  Q (•,•)  over  Q (•,•), 

c c 

matching  should  be  uniform  over  the  frequency  domain  of  interest. 

In  approximating  the  nonstationary  covariance  of  covariances 
Rc(-i»-k)  by  (C  and  substituting  in  (3.31),  the  equations  (3.23) 

are  derived.  Thus  the  covariance  method  of  linear  prediction  can  be 
derived  from  a frequency  domain  formulation  where  the  2D  spectrum  of  a 
nonstationary  signal  is  to  be  approximated  by  an  all-pole  2D  spectrum. 
Under  the  assumption  of  a stationary  signal  the  generalized  formulation 
reduces  to  the  autocorrelation  method  of  linear  prediction.  Note  that 


k 
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an  interpretation  of  the  above  spectral  error  criterion  is  not  obvious, 
due  to  the  fact  that  Qc(*,*)  and  Qc(*,*)  are  spectra  for  the  positive 

going  tail  of  a covariance  sequence  used  as  a data  sequence,  so  to  speak. 

T 

The  element  of  C C in  (3.24)  can  therefore  be  Interpreted  as  covariance 
estimates  for  a covariance  sequence.  When  one  applies  the  covariance 
method  of  linear  prediction  to  a data  sequence  of  Increasing  length,  its 
solution  approaches  the  solution  for  the  autocorrelation  method  of  linear 
prediction.  The  latter  analogy  for  covariance  sequence  prediction  cannot 
hold  because  the  covariance  sequence  estimate  for  an  ARMA(M,M~)  process 
should  be  considered  representing  a nonstationary  signal.  One  should 
therefore  be  cautious  in  extrapolating  known  linear  prediction  results 
to  the  covariance  sequence  linear  prediction  practiced  here.  We  will, 
however,  transplant  a nice  recursive  solution  procedure  for  (3.23)  to 
the  covariance  sequence  environment,  in  the  following  section. 

Assume  that  the  solution  to  (3.23)  is  found  so  that  the  poles  of 
the  ARMA(M,M“)  system  can  be  found  by  polynomial  rootfinding: 

H M 

L a z-i  - n (1-p.z-1)  ; aft  £ 1 3.34 

i-0  1 j-1  3 U 

One  way  to  solve  for  all  the  roots  of  a polynomial  is  by  the  deflation 
procedure  as  used  in  Muller's  algorithm  [CB].  An  interpretation  is  that 
the  nonlinear  system  of  equations  (3.3)  for  the  general  least  squares 
covariance  sequence  approximation  problem  is  now  concentrated  in  the 
single  algebraic  equation  of  (3.34),  which  is  a considerable  reduction 
in  complexity  indeed. 

The  covariance  sequence  approximation  procedure  proposed  here 
therefore  consists  of  the  following  steps: 


i)  Find  autoragrasslva  prediction  parameters  by  solving  the 
linear  system  o£  aquations  (3.23) 
ii)  Compute  tha  corresponding  polos  according  to  the  nonlinear 
algebraic  aquation  in  (3.34) 

ill)  Find  corresponding  linear  weights  in  a least  squares  fashion 
by  solving  the  linear  system  of  equations  In  (3.7) 

All  but  the  least  squares  fitting  aspects  of  this  procedure  have  been 
known  for  soma  two  hundred  years  in  the  function  interpolation  context 
as  Prony's  method  [PR].  Note  that  the  last  step  Is  a departure  from 
matching  conditions  used  in  the  derivation  of  the  covariance  sequence 
prediction  approach.  Note  also,  however,  the  compliance  with  the 
general  least  squares  approach,  assuming  tha  poles  are  known.  Qualita- 
tively thie  should  move  the  procedure  from  a modified  least  squares 
approach  in  the  direction  of  the  least  squares  approach.  A special  case 
arises  when  satisfies  an  ABMA(M,M~)  system  of  equations  for  then  its 
representation  is  found  exactly.  One  such  case  is  in  preserving  second 
order  statistical  properties  from  a continuous  time  system  in  a discrete 
time  Implementation.  Tha  present  method  serves  therefore  as  another 
alternative  procedure  in  tha  design  of  covariance  invariant  digital 
flltare  [PS],  [KPSS]. 

Note  that  the  system  of  equations  in  (3.7)  is  written  in  complex 
form.  Recognizing  that  the  elements  of  the  matrix  P and  vector  c are 
real  or  occur  In  complex  conjugate  pairs  enables  us  to  rewrite  the 
linear  system  of  complex  equations  (3.7)  as  a linear  system  of  real 
equations.  Complex  conjugate  columns  in  P corresponding  to  p.-pe^® 

j -18  1 

«na  Pi+1"P«  get  thereby  replaced  as  follows: 
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P A 


*2 


*L 


” 1 0 

• 

• pcosS  pslnB  • 

• • 

• 

2 2 
• p cos2B  p sin2B* 

• 

3i 

* 

A. 

>\i 

• • 

• • 

9 9 

9 9 

B 

i 

• 

9 9 

9 9 

9 9 

9 9 

9 9 

i+1 

• 

m 9 

9 9 

L L 

* p cosLB  p sinLB* 

-PrB 


3.35 


The  elements  of  Pr  and  B are  now  real.  We  furthermore  can  relate  the 
elements  of  A to  the  elements  of  B: 


Ai  " 2 5 Ai+1  " Ai  3,36 

For  real  leave  column  In  P and  element  A^  unchanged. 

3.3  An  Efficient  Recursive  Procedure  for  the  Covariance  Sequence 
Predictor 

We  will  now  take  a closer  look  at  a recently  developed  inner 
product  formulation  for  minimizing  prediction  errors  [MG].  This  pro- 
cedure will  ultimately  lead  to  fast  computational  algorithms  and  a 
sequential  test  for  stability.  The  present  development  is  based  on 
[IS],  where  methods  for  implementing  the  autocorrelation  method  of 
linear  prediction  were  presented.  The  order  m of  the  polynomial  will 
develop  recursively,  so  idf1,2, . . . ,M,  where  M is  the  preset  maximum  order 
for  the  predictor  polynomial  filter.  Let  us  define  the  prediction  error 
of  interest,  a forward  prediction  error: 


Cm(k)  ’ Ck  ”[-  1 ami  Ck-1] 


- £ a . C,  . ; a _ - 1 

ml  K-l  mO 


In  the  development  we  also  need  a backward  prediction  error: 


c~(k)  = C -[-  l b C,  , ] 

m k-m-1  1 mi  k-i 


bmi  Ck-i  5 bm,nrt-l  A 1 


The  total  forward  and  backward  prediction  errors  can  now  be  represented 


« - £ [c»00]2 

m k»M  m 


3_  **  E [c“(k)]‘ 

m k=M  m 


3.39a 


3.39b 


We  see  that  cm(k)  and  cffi(k)  can  be  considered  the  outputs  of  two 
filters  Am(z)  and  Bm(z)  having  as  common  Input  sequence  {C^} , where 
“ _i  a 

Vz)  = J;Q  V z ; am0  " 1 3-40a 


B (z)  * £ b . z 1 ; b . ■ 1 

m ml  m,m+l 


3.40b 


Let  us  now  define  the  inner  product  as  the  summation  of  the  products 
of  the  output  sequences  c*(k)  and  c~(k),  as  illustrated  in  Figure  3.1. 
This  inner  product  definition  can  be  written  as 


L m m+1 

<A_(Z)»B  (z)»  _ * £ £ a C,  £ b C. 

CTC  k-M  1-0  ml  k-i  j-1 


3.41a 
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m nri-1  L 

^ E a . E C|  . C.  . b . 
i-0  J-l  m±  k-M  1C*1  k'j  mj 


3.41b 


We  readily  recognize  the  terms  of  C C as  defined  in  (3.24),  which  c* 

be  interpreted  as  covariance  of  covariance  estimates: 

L 
E 

k=M 


1 ct-i  °k-j 


{cTc} 


ij 


* < z~±,z~J  >CTC  3.42 

Note  that  the  inner  product  represents  a number  that  depends  upon  the 
input  sequence,  the  filter  coefficients,  and  the  summation  limits  M and 
L.  The  summation  limits  determine  the  analogy  with  autocorrelation  and 
covariance  method  of  linear  prediction.  In  the  autocorrelation  method 
we  extend  the  summation  limits  indefinitely,  which  makes  it  possible  to 
interpret  the  inner  product  in  terms  of  integrals  in  the  z-plane  or  the 
frequency  domain.  Only  in  the  autocorrelation  method  is  the  problem 
of  minimizing  prediction  errors  isomorphic  with  the  problem  of  polynomial 
approximation  on  the  unit  circle  [KA].  Using  the  inner  product  defini- 
tion of  (3.41)  the  prediction  errors  of  (3.39)  can  be  expressed  as  the 
norms  of  the  corresponding  filter  polynomials: 

“m  " <Am(z)’Vz)>cTc  = I l\(z)  I I^TC  3*43a 

and 

e « <B  (z),B  (z)>  - I |B  (z) I |2_  3.43b 

m m m CTC  " m 

Minimizing  prediction  errors  can  therefore  equivalently  be  described 

as  minimizing  weighted  norm  squares  of  the  polynomials  A (z)  and  B (z) 

m m 

for  m*l,2, . . . ,M.  i?rom  the  definition  of  the  inner  product  we  see  that 
a filter  A^(z)  has  a zero  norm  if  the  output  samples  are  zero  for  kmM 
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.1 


through  k*L.  This  event  will  happen  for  finite  L if  the  input  data 
exactly  satisfies  a homogeneous  difference  equation,  whose  characteris- 
tic polynomial  is  A^(z) . The  norms  used  in  (3.43)  are  therefore  really 
pseudo  norms,  and  it  is  related  to  this  fact  that  we  cannot  give  a one 
dimensional  spectral  error  criterion  for  the  covariance  sequence  pre- 
diction error  minimization. 

If  A (z)  minimizes  a , adding  cz  j“l,2,...,m  to  the  polynomial 


m m 

must  result  In  a larger  norm  square: 


| |A  (z)  + cz  j| |2  > | |Am(z)||2T  for  j“l,2 m 

C C 


m 


T 

C C 


3.44a 


or 


2c  <A  (z) , z +c2  <z-jz^>T  Vc 

m CAC 


T 

C C 


3.44b 


If  <z  ^,z  is  non-zero,  choosing  c as 


c = “<Am(z),z  J>  /<z  \z  ■*> 

C C C C 


m 


3.45 


leads  to 


-[<Am(z),z  ]2  > 0 

C C 


m 


3.46 


If  <z  ^,z  equals  zero,  then  c = -<A  (z),z~^>  ^ leads  to  the  same 

m i 

C C 


m 


result  (3.46),  which  must  hold  for  j=l,2,...,m.  A necessary  condition 
for  minimizing  the  total  squared  prediction  error  therefore  is 

3.47 


<A _(z) , z > ~ — 0 { j=l,2, . . . ,m 

C C 


Similarly  derived  is  the  result 

0 : j — l,2,...,m 


<B  (z),z 

m cTc 


3.48 


Adding  an  orthogonal  polynomial,  say  Q(z),  shows  that  the  orthogonality 
conditions  (3.47)  and  (3.48)  respectively,  are  also  sufficient: 
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l|A.(z> + Q(z>||2  -||a(z>||2  +||q(z)||2t 

C C C C C C 

lllA_(*)||2  -a 

“ CLC  “ 


3.49 


Therefore  the  orthogonality  conditions  (3.47)  are  necessary  and  suffi- 
cient conditions  for  minimization  of  the  total  squared  forward  covariance 
sequence  prediction  errors. 

The  orthogonality  conditions  form  the  basis  for  solving  for  A (z) 

m 

and  B (z)  In  a recursive  manner.  Suppose  A (z)  and  B (z)  have  been 
m mm 

obtained.  Then  a linear  combination  of  the  form  A (z)  + k ,,B  (z)  will 

m nrt-l  in 

he  a polynomial  of  proper  order  art-1,  with  leading  coefficient  one,  and 

—1  —2  — m 

will  be  orthogonal  to  the  powers  z ,z  , ...,z  due  to  the  orthogonality 
conditions.  We  therefore  have  to  choose  k t ^ so  as  to  make  the  linear 
combination  orthogonal  to  z * , and  A^ ^(z)  will  have  been  obtained, 
and  the  forward  prediction  error  ct^ ^ will  have  been  minimized.  As 
suggested  we  will  therefore  choose 


Wz)  " Am(z)  + kart-lBm(z) 


3.50 


and  invoke  the  orthogonality  conditions  to  yield  an  expression  for 


art-1’ 


<Anrt-l(z)»z"<iafl)>  T * 0 
v C 


<A + kmtl  <Bm(z),z-<"*1>>cIc 


m 


C C 


c c 


c c 


3.51 


Since  A (z)  and  B (z)  are  both  orthogonal  to  z \...,z  m we  have 
in  m 7 

<\(*),*'tafl>>J„  ‘ <Vz)-,m(z)>„T„  ■ <1>V*)>  T 3-52 


and  also 


<»„(*). *‘<**'1)»cTc  ■ <Bm(2),Bm<z)>cT(; 


|B„(z)||  - S» 

w V 


3.53 
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Note  that  expressions  (3.52)  and  (3.53)  illustrate  that  many  inner 
product  calculations  which  might  involve  double  summations  simplify 
to  single  summation  due  to  the  orthogonality  conditions.  Substitution 


of  (3.53)  into  (3.51)  yields  the  choice  for  k ..  in  (3.50) 

IDtI 


knri-l  = ' ^ <Vz)’z"(nH’1)>  T 

mfl  3m  m CTC 


1 m T 

= - — E (C  C}  . a 
6m  i=0  mfl,i  ml 


3.54 

-1 


Recall  that  ^^(z)  is  a linear  combination  of  powers  of  z from  z to 


_ /—  j -I  \ 

z . Therefore  each  polynomial  B^(z)  is  orthogonal  to  all  lower 


order  polynomials  B^_^(z)  due  to  the  orthogonality  conditions.  The 


polynomials  (Bm(z)}  thus  form  an  orthogonal  set,  i.e. 


<B  (z),B  (z)>  _ = 0 for  m^i 

CTC 


3.55 


The  latter  property  together  with  the  orthogonality  relations  allows  a 
recursive  calculation  of  the  total  squared  forward  prediction  error  in 
each  step  of  the  recursive  process.  Start  with  Aq(z)-1  and  use  (3.50). 


Then 


A (z)  = 1 + I k.  B (z)  ; m > 0 

1-1  1 X 


3.56 


With  (3.53)  and  (3.55) 


m 


,A-<Z,-1|Ictc  ■ £ k‘ 


MAm(z)M^T  '2  <*»»!>,.  +1 


(TC 


■■  "cTc 


1 - a 


m 


3.57 


The  recursion  is  then  given  by 


a * a - - k 0 . 

m m-1  m m-1 


3.58 


| 
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Let  us  initialize  the  process  with 

-1 


Aq(z)  - 1 ; Bq(z)  - z 

Then  the  initial  total  squared  prediction  errors  are  given  by 


3.59 


a0  ” <A0*Z*»A0*Z*> _T_  " tC  CJ0,0 

w L 


Bo  " <bo(z)’bo(z)>„t„  ■ {c  C>1,1 


3.60a 

3.60b 


C C 


Now  compute 


<AQ(z),z  1>^  (3-59)  <A0(z),BQ(z)>^ 


(3.40) 


-1 


<l,z  > 


C C 


(3'-2)  <CT°>0,1  ■ {cTc>i,o 


3.61 


Substituting  the  latter  result  in  (3.54)  yields 


{cTc}io  . . {cTc}oa 


8, 


3.62 


With  the  recursive  relation  (3.50)  the  next  order  A^(z)  is  found  so 


that 


aio  • 1 ! *u  ■ *1 


3.63 


The  reduced  total  squared  prediction  error  is  found  from  (3.58): 


“1  " °0  “ klB0 


3.64 


Next  we  will  develop  recursions  for  obtaining  A (z)  for  m-l,2,...,M 


At  completion,  the  Inverse  filter  and  total  squared  prediction  error 
are  given  by 

A(z)  - A^z) 


a “ a. 


3.65 


Asstme  step  m-1  has  been  completed,  so  that  B^z)  and  8^  for  i-0,...,m~l 
are  known,  as  well  as  A (z)  and  a . To  complete  the  next  step  it  is 


J 
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necessary  to  know  Bm(z).  Recall  that  the  Bm(z)  polynomials  form  an 
orthogonal  set,  so  that  the  classical  G ram-Schmidt  orthogonallzatlon 
procedure  can  be  used  recursively: 

B (z)  - .-<■“>  - £ r 
" 1-0  ml 

Then 

-(nrt-1) 


- 0 - <*‘  - V61 


and  therefore 


3.66 


3.67 


'ml 


<z“(®fl>>B  (z)>  /B  for  B.tf) 
C C 1 


arbitrary  for  3^=0 


^ for  1*0, 1,2, . . . ,m-l  then  we  can  write 


1+1 


mi  B 


f <2'(”fl). 


3=1 


ij 


T 

C C 


i "l  ..i 


3.68 


" {C  C}nrt-l,j  bij  for  1=0 »1, •• • ,m-l 

Substituting  the  latter  in  (3.66)  yields  for  the  coefficients  of 

polynomial  B (z): 

m 


3.69 


bm,  nrt-1  1 

m-1 


m-1 


ra’J  " 'Jo  Yffll  b«  = "i-j_i  Y»*  bU 
The  total  backward  prediction  error  can  now  be  evaluated: 

».  ■ <V2>.»„<2>>  T 3-53  <r<wl) 

C C m 


3.70 


T 

C C 


3.40  -(nrt-1)  nH’1  -1 

“ <z  • 1 b«i2  T 

J"1  CTC 
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3 42  T 

- jfl  «*»,*., 


3.71 


Evaluation  o£  B (z)  and  3 are  herewith  complete  at  step  m,  and 
n in 

A^^Cz)  can  be  evaluated  from  (3.50)  and  (3.54).  The  resulting  (nrt-1) 
order  forward  predictor  coefficients  are: 

am+l,0  = 1 

amfl,i  " am,i  + knH-l  bmi  ; 


st 


amfl,m+l  ^m+1 


3.72 


The  associated  covariance  sequence  prediction  error  is  evaluated  from 
(3.58) 

°m+l  ' “m  - *Wl  em  3*73 

This  completes  a full  circle  in  the  iterative  process. 

3.4  Singularity  and  Stability  Considerations 

A Fortran  routine  implementing  the  recursive  procedures  outlined 
in  the  previous  section,  is  given  in  [MG].  These  subroutines  were  used 
in  the  simulations  performed  in  the  application  chapters  of  this  thesis. 
If  during  the  iterative  process  a norm  squared  error  equals  sero,  or 
becomes  negative,  the  iteration  can  be  terminated  since  the  result  is 
either  exceeding  the  accuracy  of  the  computer,  or  in  error.  A sero 
norm  can  only  occur  if  the  data  can  be  described  exactly  by  a linear 
combination  of  complex  exponentials,  and  this  would  Indicate  that  CTC 
in  (3.23)  is  singular.  Writing  out  (3.55)  according  to  the  definition 
gives 


m+1  n+1 

<B  (z),B  (z)>  T * l 
" " CTC  i-1  1-1 


£ b»l lc  c>l,j  b„j  ■ »•» 


In  matrix  form 


where 


BTCTC  B - DO) 


1 bll  b21 


bH-l,2 


; DO) 


Since  B is  a triangular  matrix  with  l's  on  the  diagonal  we  derive 
from  (3.75) 

T M-l 

|CAC|  - |D(1)|  - n B 3. 77 

m-0  m 

So  when  8^*0  for  some  n,  the  sequence  {C^}b  satisfies  a homogeneous 

difference  equation  of  order  n+1.  We  point  out  once  again  the  Importance 

of  the  recursive  procedure  in  that  it  detects  singularity,  or  for  that 

matter,  near  singularity,  of  the  matrix  CTC.  During  the  recursion  one 

may  get  some  idea  about  the  stability  of  the  polynomial  A (z)  as  can  be 

m 

deduced  from  (3.50): 


A (z)  - A . (z)  + k B , (z) 
m m-l  m m-l 


tii 


The  coefficient  with  the  highest  power  of  z-1  yields: 


11  H>  i> 
l-l  ml 


Therefore,  if  |k  | > 1,  then  at  least  one  of  the  poles  {p  }®  has 

radius  larger  than  one,  and  the  polynomial  A (z)  represents  an  unstable 

n 
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system.  Alternatively,  we  can  state  that  a necessary  condition  for 
the  stability  of  Am(z)  is: 

l*J  < 1 3.80 

The  coefficients  k above  relate  to  both  polynomials  A (z)  and  B (*) 

m hi  hi' 

however,  and  are  therefore  not  sufficient  to  determine  A (z).  One 

m 

can  find  a similar  but  uniquely  related  set  of  coefficients  k when  the 

n 

polynomials  Bm(z)  are  related  to  the  polynomials  Affl(z)  by 

Vz)  “ z'(m+1)  Am(l/z)  3.81 

The  latter  equality  is  valid  when  {CTC Is  a function  of  (i-j)  only,  as 

in  the  autocorrelation  method  of  linear  prediction,  under  the  assumption 

of  stationarity . In  making  a similar  definition  however  we  can  derive 

a set  of  reflection  coefficients  {1^}“  that  uniquely  relates  to  the 

polynomial  A (z).  Let  us  define 
m > 

(^(z)  = z-(mfl)  Am(l/z)  i Bffl(z)  3.82 

Combination  of  (3.82)  and  (3.78)  yields 

Qm(z)  * z"1[,cmAm-l(z)  + Vl(z)]  3.83a 

or 

Vl(z)  " zQm(z)  - Km  Am-l(z)  3-838 

Substituting  in  (3.78)  and  solving  for  A ,(z)  in  terms  of  A (z) 

m— l in' 

gives : 


Am-l(z) 


Am(z)-ZKmVz) 


1-K 


+ 1 


With  (3.82) 


3.84a 


Vl(z) 


Am(z)-2~mKn)Am(1/z> 

1-K2 

m 


5 lKml  * 1 


3.84b 
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or 


m~l 

l 3.  ■ , z 

i=0  m_1»i 


m -i  m - i 

Ea.z-ic  E a ,z 

-i  i=0  mi  m i-0  m’m_i 


1-k 


l<J  * 1 3.84c 


m 


From  the  higher  order  polynomial  A (z)  and  k we  can  therefore 


m 


m 


recursively  find  A , (z): 

m-i 


a -k  a 


mi  m m,m-i  , , . 

m-l  i ™ 2 * i“0,l...,m-l  , 

* 1-k 

m 


< I4  1 
m 


3.85 


and  thereby  k 


m-i ' 

Km-1  am-l,m-l 


3.86 


The  latter  procedure  provides  an  efficient  recursive  procedure  to  test 

the  magnitude  of  the  reflection  coefficients  Km>  corresponding  uniquely 

to  the  filter  polynomial  A^Cz),  at  every  step  of  the  recursive  process 

to  find  the  optimal  linear  covariance  sequence  predictor.  A necessary 

and  sufficient  condition  for  stability  of  the  polynomial  Affl(z)  is  [MG] 

| K | < 1 V m 3.87 

m 

In  this  section  a recursive  procedure  has  been  given  to  solve  for 
the  best  linear  covariance  sequence  predictor  (3.23).  The  order  of 
the  predictor  polynomial  is  increased  in  every  step.  Efficient  proce- 
dures have  been  given  to  evaluate  the  stability  of  the  corresponding 
inverse  filter  as  well  as  the  residual  total  squared  covariance  sequence 
prediction  error  for  that  step.  One  can  therefore  stop  increasing  the 

order  when  Instability  is  indicated  and  say  that  A .(z)  is  the  best 

m— l 

linear  predictor  of  maximum  order  that  has  a stable  Inverse.  Also, 

whenever  6_  is  numt-rically  equal  to  zero  we  have  arrived  at  the  best 

m 

linear  predictor  of  minimum  order.  The  latter  property  enables  us  to 
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find  the  desired  solution  even  if  C C in  (3.23)  is  singular.  Note  that 

we  assumed  the  maximum  order  of  the  underlying  system  to  be  M.  Such  an 

assumption  requires  an  educated  guess  or  a restriction  due  to  implements- 

tion  requirements.  In  the  above  procedures  however  it  is  not  necessarily 

a disadvantage  to  choose  the  maximum  order  somewhat  too  high,  as  during 

the  recursive  process  the  residual  prediction  error  a provides  an 

m 

indication  of  the  resulting  performance.  The  behavior  of  {a^}  thus  pro- 
vides a built-in  estimate  of  the  appropriate  system  order. 

Similarities  of  the  present  pole  identifier  with  the  covariance 
method  of  linear  prediction  of  speech  are  clear.  In  the  latter  applica- 
tion the  need  arises  for  formant  estimation.  A damped  sinusoidal 
component  of  the  vocal  tract  acoustic  impulse  response  is  a formant. 

The  analogy  with  the  damped  exponentials  {p^}  in  our  covariance  sequence 
model  of  (3.1)  is  obvious.  Formant  estimation  for  speech  signals  is 
therefore  an  analog  for  pole  identification  or  estimation  for  covariance 
sequences.  The  accuracy  of  formant  estimation  obtained  with  the  covari- 
ance method  of  linear  prediction  has  turned  out  to  be  vastly  superior 
to  the  results  derived  with  the  more  popular  autocorrelation  method  of 
linear  prediction  [CL].  In  particular  when  only  short  data  sequences 
were  available  this  supremacy  is  pronounced.  The  desirability  of  the 
linear  prediction  approach  for  formant  estimation  is  qualitatively 
justified  by  the  fact  that  in  signal  representation,  the  most  efficient 
basis  set  (or  set  of  approximating  functions)  is  the  one  which  most 
closely  matches,  in  some  sense,  the  desired  characteristics  of  the 
signal.  The  latter  statement  reflects  remarks  made  at  the  closing  of 
Section  5.2.  We  indicated  there  that  the  proposed  covariance  sequence 
model  will  fit  certain  processes  with  the  minimum  nun&er  of  coefficients 
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because  of  its  underlying  structure.  In  the  exposition  for  the  design 
of  a certain  class  of  recursive  digital  filters  this  philosophy  will 
play  an  important  role.  See  also  [GRT]. 

I ' I 
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4 COVARIANCE  SEQUENCE  APPROXIMATION  FOR 


PARAMETRIC  SPECTRUM  ANALYSIS 


An  introduction  is  given  in  Section  4.1  that  concentrates  on 
parametric  spectrum  analysis  and  why  it  might  be  beneficial  to  take 
such  an  approach.  The  section  after  contains  a derivation  of  the  para- 
metric spectral  density  estimate  associated  with  the  covariance  sequence 
approximation  procedure  of  Chapter  3.  Its  properties  and  peculiarities 
are  examined  as  well.  Section  4.3  extends  the  covariance  sequence 
approximant  and  the  associated  spectral  density  estimate  so  as  to  include 
the  analysis  of  harmonic  processes  in  additive  noise.  Sections  4.4  and 
4.5  give  results  for  the  covariance  sequence  approximant  applied  to 
known  and  estimated  covariance  sequences,  respectively.  Application  of 
the  basic  idea  to  the  fitting  of  a known  but  nonrational  Gaussian  spec- 
tral density  indicates  the  technique  is  robust. 

4.1  Parametric  Spectrum  Analysis  Introduction 

In  order  to  see  parametric  spectrum  analysis  procedures  in  the 
proper  perspective  we  will  briefly  review  some  ideas  pertaining  also  to 
nonpar ame trie,  spectrum  analysis.  Let  us  first  note  that  our  interest 
will  be  in  spectral  estimation  for  sampled  data.  For  estimation  of 
spectra  from  continuous  time  systems  this  means  that  the  resulting  spec- 
tral estimate  is  an  infinite  sum  of  translates  of  the  continuous  time 
spectrum.  The  socalled  aliasing  effect  will  be  ignorable  for  a proper 
choice  of  the  sampling  frequency.  Exact  determination  of  either  covari- 
ance function  or  spectrum  would  require  a perfectly  measured,  infinitely 
long  realization  and  infinite  wordlength  computation  for  all  wide  sense 
stationary  processes  with  nondeterminlstic  components.  Both  require- 
ments are  impractical.  Approximate  determination  however  immediately 
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brings  up  the  questions  of  how  long  the  observation  interval  should  be, 
which  computational  approach  should  be  used  and  finally  how  much 
confidence  one  can  have  in  the  resulting  estimate.  The  answers  are 
relatively  simple  for  nonparametric  spectrum  estimation  [BT],  but  unsatis- 
factory because  of  the  required  observation  length  for  highly  precise 
estimates.  The  finite  observation  interval  has  the  effect  that  the 
resulting  spectral  estimate  is  at  best  a convolution  of  the  continuous 
time  spectrum  with  the  l.anczos  spectral  window  corresponding  to  a finite 
length  observation  sequence.  The  latter  has  the  effect  of  reducing  the 
spectral  resolution.  It  is  imperative  therefore  that  one  realizes  that 
a finite  length  data  sequence  yields  no  more  than  estimates  for  the 
covariance  function  and  the  spectrum,  and  that  consequently  these  esti- 
mates exhibit  bias  and  variance.  In  parametric  approaches  however 
observation  window  and  spectrum  estimate  are  deconvolved  in  some  sense 
by  implicit  extension  of  the  data  sequence. 

Many  methods  for  spectrum  estimation  have  been  proposed  over  the 
years,  both  parametric  and  nonparametric  ones.  Methods  based  on  finding 
a spectrum  estimate  for  the  observed  data  directly  have  received  an 
impulse  due  to  the  computational  advantages  of  the  FFT  algorithm.  The 
underlying  assumption  of  periodicity  was  far  outweighed  by  the  conve- 
niences of  the  method.  A finite  data  sequence  yields  a finite  covariance 
sequence  and  the  corresponding  spectrum  has  finite  spectral  resolution 
due  to  the  finite  observation  window.  The  resolution  properties  of  the 
conventional  spectral  analysis  methods  are  chiefly  due  to  the  assumptions 
that  are  made  concerning  the  data  outside  of  the  observation  interval. 

The  desire  to  be  consistent  with  the  available  covariance  seouenee  and 
at  the  same  time  he  noncommittal  about  the  unobserved  data,  has  led 
recently  to  the  maximum  entropy  spectral  estimator  |B1],  also  known  as 
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a high  resolution  spectral  estimator.  The  maximum  entropy  estimator  can 
be  derived  from  maximum  entropy  considerations  for  Gaussian  processes, 
and  it  was  shown  later  that  the  corresponding  extrapolation  of  the 
covariance  sequence  is  equivalent  to  least  squares  fitting  of  an  auto- 
regressive model  to  the  process  [VPB] . The  high  resolution  is  achieved 
by  parametric  modeling  and  implicity  extending  the  available  data.  The 
Burg  algorithm  does  not  explicitly  make  an  estimate  of  the  covariance 
sequence  lags.  Instead  it  operates  directly  on  the  observed  data 
samples  and  minimizes  the  sum  of  forward  and  backward  prediction  errors 
simultaneously  [B2].  It  has  been  found  that  this  procedure  is  quite 
sensitive  to  the  initial  phase  of  observed  sinusoids  [CS].  A least 
squares  estimator  minimizing  the  same  error  criterion  as  in  the  Burg 
algorithm,  but  based  on  explicit  computation  of  the  corresponding  covari- 
ance sequence,  was  found  to  be  much  more  stable  with  respect  to  the 
initial  phase  of  the  signal  [UC] . The  high  resolution  spectral  estimators 
have  been  very  successful  in  resolving  pronounced  spectral  peaks,  but  a 
word  of  caution  is  not  out  of  place.  The  notion  of  high  resolution  would 
seem  to  indicate  that  all  spectral  features  can  be  resolved  well.  It  is 
known  however  that  influential  zeros  in  a spectrum  are  not  modeled  very 
well  by  the  underlying  all  pole  model  of  any  autoregressive  spectral 
estimator.  We  rust  therefore  recognize  that  the  term  high  resolution 
spectral  estimation  has  chiefly  pertained  to  the  resolution  of  spectral 
peaks.  The  smooth  nature  of  the  corresponding  spectrum  estimate  is  due 
to  the  fact  that  it  is  the  spectrum  of  a process,  for  which  the  available 
data  merely  represents  a realization.  Such  a process  spectrum  will  of 
course  look  much  smoother  than  a spectrum  estimate  based  on  the  realiza- 
tion directly. 
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In  the  parametric  spectrum  estimation  approach  of  this  thesis  we 
pose  the  basic  question  of  parametric  spectrum  analysis:  Which  system, 
driven  by  white  noise,  is  likely  to  have  generated  the  present  finite 
realization?  The  spectrum  associated  with  that  system  is  called  the 
spectral  estimate.  For  many  applications  one  furthermore  wishes  to  have 
a low  parameter  representation  for  purposes  of  bandwidth  efficient  com- 
munication as  in  linear  predictive  coding  of  speech  [MG] . Maximum 
entropy  spectral  analysis  and  other  autoregressive  spectrum  estimators, 
lead  to  a spectral  density  ploL  and  an  AR  model,  but  do  not  readily  give 
frequency  and  corresponding  power  estimates.  Reasonable  approximations 
have  been  suggested  for  low  noise  cases  with  known  covariances  [ JA] . 

The  latter  method  based  on  approximating  complex  residues  yielded 
improved  estimates  over  those  based  on  integration  of  the  spectral 
density  around  peaks  [LA]. 

One  of  the  most  serious  problems  in  parametric  spectrum  analysis  is 
the  determination  of  an  appropriate  order  for  the  underlying  model.  If 
in  the  autoregressive  spectrum  estimators  too  few  parameters  are  used  a 
highly  smoothed  spectral  estimate  results.  If  on  the  other  hand  too 
many  parameters  are  determined,  spurious  peaks  will  turn  up.  This 
order  determination  problem  is  most  pronounced  for  estimation  from  short 
realizations.  In  the  proposed  approach  of  this  thesis  spurious  poles  ca. 
turn  up  in  the  initial  stage  of  the  approximation,  but  the  corresponding 
linear  combination  weights  provide  a correction  mechanism  to  assess  their 
proper  importance.  For  maximum  entropy  and  other  autoregressive  spectral 
estimators  order  determination  criteria  have  been  proposed  and  in  many 
instances  used  quite  successfully  [AK] , [PA].  In  the  proposed  spectral 
estimator  the  order  is  increased  recursively  and  a simultaneous 


indication  of  the  covariance  sequence  prediction  error  points  out  if  the 
order  is  adequate  before  reaching  a preset  may -I  mum  order.  An  error  cri- 
terion as  Akaike's  starts  to  increase  at  some  order  and  intuitively 
reflects  the  danger  of  using  an  ever  increasing  model  order  due  to  the 
fact  that  with  high  enough  order  even  extremely  noisy  data  can  be  fitted 
exactly.  This  philosophy  led  in  the  system  identification  and  time 
series  modeling  area  to  the  Croup  Method  of  Data  Handling  [IV].  In  the 
GMDH  method  model  parameters  are  based  on  one  data  set  and  model 
appropriateness  is  judged  from  another,  preferably  disjoint  data  set. 

The  assumption  of  sinusoids  in  white  noise  led  to  the  covariance 
sequence  model  associated  with  the  Pisarenko  decomposition  [PI]. 
Pisarenko  points  out  that  the  maximum  entropy  spectrum  for  harmonic 
processes  in  additive  white  noise  is  a smoothed  version  of  the  Pisarenko 
decomposition  result.  As  no  autoregressive  spectral  estimator  has  any 
zeros  in  its  spectrum,  it  is  obviously  the  proximity  of  the  zeros  to  the 
unit  circle  that  determines  the  resolution  power  of  spectral  estimators 
for  harmonic  processes.  The  latter  property  indicates  the  usefulness 
of  ABMA  spectral  estimators,  several  of  which  have  been  proposed  [G2], 
[KPSS],  [K],  [TS].  The  minimum  residual  criterion  in  the  latter 
determines  the  zeros  first  with  an  iterative  Powell  minimization  tech- 
nique and  thereafter  performs  a further  whitening  procedure  on  the 
filtered  signal  resulting  in  a linear  system  of  equations  for  finding 
the  poles.  When  the  actual  signal  spectrum  has  Influential  zeros,  quite 
a few  poles  are  needed  in  an  all  pole  model  to  achieve  a good  approxima- 
tion. The  introduction  of  zeros  will  then  certainly  allow  a more 
efficient  representation. 
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Part  of  the  significance  of  ARMA  modeling  of  course  lies  in  the 
fact  that  sampling  of  a continuous  AR  process  of  order  M results  in  a 
discrete  ARMA(M.H-l)  [BA].  That  property  holds  for  any  strictly  proper 
continuous  time  process  of  order  M [P ] . Furthermore,  any  continuous 
spectral  density  function  can  be  approximated  arbitrarily  closely  by  an 
ARMA(M,N)  process,  by  proper  choice  of  the  parameters.  While  this  is 
also  true  of  spectral  densities  of  finite  MA  and  finite  AR  processes, 
a model  with  a rational  spectral  density  will  almost  Invariably  lead  to 
a better  fit  with  fewer  parameters  [K] . The  latter  property  has  been 
referred  to  as  one  of  parsimony  [BJ]. 

4.2  Parametric  Spectrum  Estimation  Based  on  Covariance  Sequence 
Approximation 

Recall  that  based  on  a finite  covariance  sequence  sample  { 
we  find  an  approximating  sequence 

\ - X p±lk|  *k  4-1 

where  {A^.p^}  are  real  or  occur  in  complex  conjugate  pairs.  Since  the 
spectrum  is  uniquely  related  to  the  covariance  sequence,  we  may  find 
the  approximating  spectrum  as  follows 

SU-e^8)  = HR^  = l \ z~k 
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The  region  of  absolute  convergence  is  the  annulus  determined  by  the 
radii  of  max  | | and  its  inverse.  The  approximating  spectrum  does  not 


exist  as  a uniform  limit  If  any  of  the  complex  exponentials  has  radius 
greater  than  or  equal  to  one.  This  situation  would  correspond  to  a non- 
minimum  phase  covariance  sequence  prediction  filter.  As  explained  in  the 
previous  chapter  this  can  be  prevented  from  happening  in  the  recursive 
determination  of  the  filter.  We  further  note  here  that  the  convergence 
of  S(z)  has  to  be  examined  separately  if  the  maximum  pole  radius  equals 


Let  us  examine  some  properties  of  this  spectrum  estimate.  From 
(4.2)  it  can  be  written  as  follows 


-1 

n (z-p  >(z  -p  ) 
i-1  1 1 


A(z)A(z_1) 


where  A(z)  is  the  denominator  polynomial  of  the  underlying  ARMA(M,M~) 
system.  The  numerator  N(z)  is  given  by: 

M M 1 

N(z)  " E Vl'Pi) <1+V  n <z-Pk><*  -pk>  4.3b 

-l  ^ 

Obviously  N(z)  = N(z  ),  so  if  z is  a root  of  N(z),  then  1/z  is  also 

O o 

a root.  Let  us  define  the  following  polynomial  [K] 

<Kz)  - z(M_1)N(z)  4.4 

this  polynomial  is  of  degree  2(M-1)  in  z.  Again  if  z is  a root  then  1/z 

o o 

is  a root,  and  the  polynominal  can  therefore  be  written 

M— 1 

<Kz)  “ Const  n (z— Zj ) (zZj— 1)  4.5 


From  (4.4) 


N(z)  = z (z) 

= Const  M~1  (z-z^)(z_1-Zj) 

= B(z)B(z  l) 
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B(z)  contains  one  root  of  each  of  the  pairs  (z  , 1/z  ) and  the  Const. 

J J 

depends  on  the  choice  of  roots.  For  a given  polynominal  N(z)  we  can 
find  several  polynominals  B(z),  for  each  of  which  the  approximate 
spectrum  can  be  written  in  factored  form: 


A(z)A(z_1) 


1 


As  a result  we  say  that  there  exist  several  ARMA(M,M~)  systems  with 
transfer  function 

HARMA(z)  = A(i}  4,8 

A 

and  spectrum  S(z)  provided  B(z)  has  real  coefficients.  The  denominator 
polynomial  has  as  its  coefficients  the  predictor  polynominal  coefficients 
identified  in  the  previous  chapter  and  they  are  all  real.  Also  from 
(4.3b)  it  is  clear  that  N(z)  in  a mirror  image  polynomial  with  real 
coefficients  due  to  the  fact  that  {A^,p^}  are  real  or  occur  in  complex 
conjugate  pairs.  For  N(z)  we  then  have  the  following  implications 

N(z  ) =*  0 -*•  N(z*)  = 0 4.9a 

o o 

N(zq)  = 0 -►  N(zo_1)  = 0 4.9b 

J9 

A zero  inside  the  unit  circle,  say  zq  = pe  ; 0 < p < 1,  will  therefore 
mean  that  the  following  four  zeros  exist 

je  -je  , -je  , je 

{pe  °,  pe  ° , - e ° , ^ e °}  4.10 

We  have  seen  that  B(z)  in  (4.6)  contained  one  root  of  each  of  the  pairs 

je  , je 

(Zj , 1/Zj),  and  we  can  easily  see  now  that  choosing  pe  ° and  — e 0 
will  not  lead  to  an  AKMA  numerator  polynomial  B(z)  with  real  coefficients. 
Therefore,  if  we  desire  an  ARMA  system  with  real  coefficients  it  is 
imperative  to  choose  the  zeros  in  complex  conjugate  pairs.  Of  the  zeros 
in  (4.10)  only  one  free  choice  can  be  made.  If  a zero  of  N(z)  occurs 


inside  the  unit  circle  and  is  real,  another  zero  will  occur  outside  the 


unit  circle,  with  inverse  radius.  The  necessary  conditions  (4.9)  are 
satisfied  even  when  only  one  extra  zero  occurs.  Either  one  of  these 
real  roots  can  be  chosen  for  the  numerator  polynomial  B(z). 

Let  us  now  investigate  the  occurrence  of  a non-repeated  zero  on 

je 

the  unit  circle:  zq  = e . In  total  compliance  with  (4.9)  we  have 
as  the  only  other  zero  z*  = e ° . Again  either  one  of  these  roots 

can  be  chosen  for  the  polynomial  B(z).  This  complex  zero  does  not  occur 
in  a complex  conjugate  pair  however,  and  the  resulting  ARMA  system  has 
complex  numerator  coefficients.  For  B(z)  to  have  real  coefficients  the 
zeros  of  N(z)  on  the  unit  circle  must  occur  in  multiples  of  two.  The 
latter  requirement  turns  up  invariably  in  spectral  factorization  problems 
[VT]-  Violation  of  this  requirement  indicates  that  the  corresponding 
spectrum  estimate  is  negative  for  some  frequencies.  The  conditions  of 
(4.9)  which  are  satisfied  in  the  spectrum  estimate  of  (4.2),  cannot 
guarantee  the  spectrum  estimate  to  be  nonnegative.  Equivalently,  the 
covariance  sequence  approximant  is  not  necessarily  nonnegative  definite, 
in  which  case  it  is  not  a covariance  sequence.  It  should  be  noted  that 
negative  spectrum  estimates  are  not  uncommon,  and  may  occur  in  particular 
when  the  observations  do  not  exactly  satisfy  the  underlying  model.  One 
can  then  take  the  spectral  estimate  to  be  zero  [Gl],  or  alternatively 
plot  absolute  value  of  S.  If  the  interest  is  purely  in  a spectral 
estimate  defining  the  estimate  to  be  zero  where  it  is  negative  seems  at 
least  a reasonable  thing  to  do.  If  there  is  interest  in  the  complex 
ARMA(M,M  ) associated  with  the  spectral  estimate,  then  taking  the 
absolute  value  is  more  appropriate. 

Since  the  denominator  polynomial  A(z)  has  real  coefficients  one  can 
find  the  corresponding  complex  linear  weights  A^  in  the  model  of  (4.1) 
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i. 


under  the  condition  that  the  resulting  covariance  approximant  be 
nonnegative  definite  or  equivalently  the  spectrum  estimate  be  nonnega- 
tive. It  is  fairly  easy  to  give  sufficient  conditions,  such  as 

A±  real,  >0  ; 1 - 1 H A. 11 

Even  a substantially  less  restrictive  condition  can  be  given  by  direct 

application  of  a continuous  time  restriction  [Y] : 

a 

|<l> . j 1 arctan  -i  4.12 

Pi 

-j*i  it  ' 

Ai  “ rie  ; Pi  " pie 

“i  " * *n  P± 

The  above  restrictions  would  assure  the  nonnegative  definiteness  of 
each  of  the  individual  covariance  sequences  corresponding  to  a complex 
conjugate  pole  pair.  It  is  very  difficult  to  derive  necessary  and 
sufficient  conditions  for  the  general  model  of  (4.1)  to  result  in  a 
nonnegative  spectrum  estimate.  Imposition  of  these  conditions  would 
immediately  result  in  a nonlinear  programming  problem,  which  is  what  we 
are  trying  to  avoid. 

4.3  Harmonic  Processes  in  Additive  Noise 

An  important  class  of  processes  is  formed  by  a linear  combination 
of  random  amplitude,  random  phase  sinusoids  in  additive  noise.  In  case 
the  additive  noise  is  white  this  is  exactly  the  model  that  underlies 
Pisarenko's  method  of  spectral  decomposition  [PI).  Strictly  speaking 
one  cannot  define  the  spectrum  for  a random  amplitude,  random  phase 
sinusoid.  Consequently,  unless  one  proceeds  with  some  care,  the  spectrum 
of  a random  amplitude,  random  phase  sinusoid  has  no  well-defined  meaning. 
Therefore,  let  us  proceed  as  follows.  Define  the  process 

\ - asin(k61  + <fr)  ; k-0,  +1,...  4.13 


A 
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with  a:  N(0,a  ),  say,  and  U[0,2ir].  The  covariance  of 


* a coske^ 


k - 0,  ±1,.., 


This  covariance  sequence  is  neither  l nor  i , so  for  the  time  being 
the  formal  spectrum  with  double  poles  at  z = e , 


2 j6l  j0l 

x) 

te-*°-e3  Ste-J6-. 


2 -J01  -J0l 

2l  (1-e  1)(l+e 

2 , je  -J0iw  -j0  -J0i, 

(e  -e  ) (e  J -e  ) 


has  only  formal  meaning.  To  give  it  a well  defined  sense,  and  to  show 
how  it  may  be  obtained  as  a limiting  case  of  the  spectrum  given  in  (4.2), 
define  the  covariance  sequence 


■ P^l  o^coskG^  . 


0 < p < 1 

±J®i 


The  corresponding  spectrum,  with  poles  at  z = pe  , is  well  defined. 


From  (4.2) 


S (z=eJ0)  = ^ 

P 2 


I j°l  j0l 

2 (1-pe  ) (1+pe  X) 


-J®, 

(1-pe  i)(l+pe  L) 


Furthermore,  Sp(z=eJ  ) may  be  written  [LLS] , [F] 


II 

sp(z=ej0)  = | Pp  (9“$)S^(z*e^)d$ 


Here  P (0)  is  the  Poisson  kernel: 
P 

, 2 

p (C)  = - 

2ir(l-2pcos0  + p ) 


4.18a 


4.18b 


By  the  delta  function  character  of  the  Poisson  kernel  sequence , we  have 
sp(z-eJ°)  converging  to  S^r.-e^0)  at  all  points  of  continuity  fl.]: 

S(j(*'c-,°)  ->  S)(/.-u-,())  4.19 

pH 
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Thus  the  formal  spectral  density  S^*)  for  a random  amplitude,  random 
phase  sinusoid  may  be  obtained  as  a uniform  limit  (as  ptl)  of  the  spectra 
defined  by  our  covariance  sequence  approxlmant.  The  above  procedure 
generalizes  to  a linear  combination  of  random  amplitude,  random  phase 
sinusoids  in  an  obvious  way,  by  the  independence  of  the  amplitude  and 
phase  variables.  A m4x  of  discrete  and  continuous  spectral  density 
functions  can  now  be  represented  by  the  covariance  sequence  approxlmant 
of  (4.1)  by  virtue  of  the  independence  of  the  discrete  and  continuous 
components  as  in  the  Wold  decomposition  theorem. 

Due  to  symmetry  properties  we  are  mainly  concerned  with  the  positive 
going  part  of  the  covariance  sequence.  The  covariance  sequence  approxi- 
mant  for  nonnegative  lags  is  easily  recognized  as  the  general  solution 
of  a linear  homogeneous  difference  equation  of  order  M with  real 
coefficients  and  a characteristic  equation  with  simple  roots.  The  first 
M lag  values  are  determined  by  the  initial  conditions  and  for  lag  values 
M and  up  the  behavior  of  the  solution  is  determined  by  the  characteris- 
tic roots.  Viewing  the  covariance  sequence  approxlmant  as  the  solution 
of  a linear  homogeneous  difference  equation  we  readily  recognize  its 
validity  for  all  characteristic  roots,  whether  in,  on  or  outside  of  the 
unit  circle. 

In  order  to  show  how  white  noise  can  be  represented  by  the  covariance 
sequence  approxlmant  of  (4.1)  let  us  assume  a first  order  autoregressive 
process  driven  by  white  noise  e:  (0,1).  Assume  further  that  the 
autoregressive  coefficient  is  positive  real.  The  resulting  noise  process 
nk  is  described  by 

"k'1  Vi*ct  4-20 


70 


Multiplying  by  o and  taking  expected  values,  a recursive  relation 

for  the  noise  covariance  sequence  is  found 


c^-pcj  -0 


V±  > 0 


Consequently  we  can  find  the  positive  going  solution  to  (4.21)  and 

extend  it  for  negative  lag  values  by  the  symmetry  property  to  yield  the 

explicit  expression  for  the  noise  covariance  sequence: 

2 

“ — ^~2  Vk  4.22 

1-P 


Note  that  if  we  choose  0 < p < 1 then  the  sequence  of  sequences  {C?}  , 

indexed  by  p,  converges  in  to  the  dirac  delta  sequence  {o£6(k)}  ^ , 
as  p approaches  zero.  This  uniform  convergence  furthermore  implies  that 
the  spectrum  Sp(0)  associated  with  (4.22)  converges  uniformly  and 
absolutely  to  1 for  0e  [— it ,ir  ] . With  the  radius  of  the  noisepole  p 
approaching  zero  we  therefore  have  an  excellent  approximation  to  a white 
noise  process.  As  the  covariance  sequence  in  (4.22)  is  seen  to  be  a 
special  case  of  the  general  covariance  sequence  approximant  of  (4.1)  with 
positive  real,  we  can  state  that  the  covariance  sequence  for  the  white 
noise  process  is  a limiting  case  of  the  covariance  sequence  approximant. 
Applying  the  limiting  arguments  for  sinusoids  and  white  noise  simultane- 
ously, Pisarenko's  decomposition  technique  becomes  a special  case  of  the 
covariance  sequence  approximant  proposed  in  Chapter  2 of  this  thesis.  We 
furthermore  note  that  both  limiting  arguments  above  resulted  in  positive 
real  linear  combinations  of  complex  exponentials  so  that  sinusoids  in 
white  noise  also  form  a limiting  case  for  the  generalized  wide  sense 
stationary  process  of  Section  2.5. 

With  two  complex  poles  per  sinusoid  the  covariance  sequence 
approximant  implicitly  models  a sinusoid  in  white  noise  with  an 
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ARMA(3,2)  model.  In  view  of  the  ARMA  model  for  harmonic  processes  in 
white  noise  derived  by  Ulrych  and  Clayton  [UC],  the  present  modeling 
approach  is  not  the  most  parsimonious  one.  We  note,  however, that  in  the 
present  approach  no  eigenproblem  needs  to  be  solved.  Furthermore,  in 
the  most  parsimonious  ARMA(p,p)  model,  numerator  and  denominator  poly- 
nomial are  equal,  due  to  the  fact  that  the  model  describes  the  process 
as  if  the  additive  noise  were  its  input.  The  only  information  in  the 
ARMA(p,p)  model  therefore  concerns  the  pole  locations  exclusively.  We 
emphasize  again  that  sinusoids  in  white  noise  represent  a special  case 
for  the  covariance  sequence  approximant.  The  sinusoids  happen  to  corre- 
spond to  poles  on  the  unit  circle  and  the  white  noise  corresponds  to  an 
autoregressive  pole  at  zero.  In  the  sections  to  follow  we  will  examine 
the  estimator  results  for  known  covariance  sequences  first,  and  there- 
after for  estimated  covariance  sequences.  We  will  furthermore  experi- 
mentally test  the  robustness  of  the  covariance  sequence  approximant 
against  nonrational  spectra  by  applying  it  to  the  known  covariance 
sequence  for  the  Gaussian  spectral  density. 

4.4  Numerical  Results  for  Known  Covariance  Sequences 

In  this  section  we  examine  the  performance  of  the  covariance 
sequence  approximant,  and  the  corresponding  spectral  density  estimate, 
when  the  given  finite  length  covariance  sequence  samples  are  known.  The 
resulting  covariance  sequence  approximant  is  given  in  numerical  form  by 
real  and  imaginary  parts  of  the  linear  combination  factors  A^  and  by 
radius  and  angle  (In  radians)  for  L lie  complex  exponentials  p^.  Corre- 
sponding to  the  covariance  sequence  approximant  is  the  parametric  spec- 
tral density  plot  of  10  log^l S(z-e^9)  | ; 0 ± 0 <_  1 (in  n radians). 

Recall  that  absolute  value  is  necessary  due  to  the  fact  that  our 
"spectral  density"  estimate  may  be  negative  for  some  frequencies.  For 
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poles  with  radius  within  1^^  of  1 the  discrete  spectral  density  Is 
superimposed  on  the  spectral  density  plot  and  depicted  as  a solid  verti- 
cal line  at  the  frequency  0^  determined  for  the  sinusoid.  The  maximum 
value  of  the  spectral  line  is  given  by  10  log^g2|AjJ  , which  indicates 
the  power  of  the  sinusoid  when  is  real.  In  every  example  a maximum 
order  M was  chosen  a priori.  Unless  otherwise  indicated  poles  are 
determined  from  lags  M and  up,  whereas  linear  combination  weights  are 
determined  from  all  lags  given.  All  computations  were  performed  in 
single  precision  arithmetic  on  the  CYBER  171  at  Colorado  State  University. 
Example  1 A single  sinusoid 

Given  C = cosTr/4.k  for  lags  k = 0,...,4 

K- 

Maximum  predictor  order:  M = 2 
Covariance  Sequence  Approximant: 

R^  = (.5000  + j -488510  _14)  (1.000  exp  j .7854)^  + c.c. 

where  c.c.  in  the  sequel  denotes  complex  conjugate  of  the 
term  in  front  of  it. 

The  identified  frequency  is  correct  to  four  decimal  places.  Note 
also  that  the  power  estimate  given  by  Rq  is  correct  to  at  least  four 
decimal  places.  The  imaginary  part  of  the  linear  weights  is  truly 
negligible  as  compared  to  the  real  part,  which  is  to  be  expected  in  view 
of  the  representation  in  (4.19).  See  Figure  4.1  for  the  spectral  plot. 

Example  2 A single  sinusoid  in  white  noise 

A.  Given  C^  * cosn/4.k  + .0001  6(k)  for  lags  k = 0,...,6  (SNR=+40dB) 
Maximum  predictor  order:  M = 3 
Covariance  Sequence  Approximant: 

Rk  = (.5000  - j .563110"10) (1.000  exp  J .7854) ^ + c.c. 

+ (.100010-3  )(.568410~5  exp  -j  3.142) I k ^ 


0.60 


0.80 


1.00 


FREQUENCY  (in  tt  rad) 


Figure  4.1  Spectral  estimate  for  sine  from  5 known  covariance  lags. 
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The  identified  frequency  and  power  are  correct  to  at  least  four 
decimal  places.  Note  here  that  the  noisepole  results  in  an  excellent 
first  order  autoregressive  noise  approximation  to  the  given  white  noise. 

Note  also  that  the  white  noise  power  has  been  determined  to  at  least  four 
decimal  places  accurately.  See  Figure  4.2A  for  the  corresponding 
spectral  plot. 

B.  Given  = cosir/4.k  + 1000. 6 (k)  for  lags  k =i  0,...,6  (SNR*-30dB) 

Maximum  predictor  order:  M = 3 

] 

Covariance  Sequence  Approximant : 

-7  Ikl 

= (.5000  + j -17681()  )(1.000  exp  j . 7854) 1 1 + c.c. 

+ (.iooo1Q4  H.iooo10"9  )lkl 

The  identified  frequency  and  power  are  correct  to  at  least  four 
decimal  places.  We  nake  an  important  computational  comment  here.  With 
the  occurrence  of  a noisepole  with  radius  less  than  10  the  final 
autoregressive  coefficient  is  very  small.  To  avoid  convergence  problems 
in  the  root  finding  algorithm,  which  computes  z * roots,  the  noisepole 
is  arbitrarily  set  to  pn  = I.^q  ^ whenever  the  final  autoregressive 
coefficient  is  smaller  than  in  magnitude.  Note  the  accurate 

determination  of  the  noisepower.  See  Figure  4.2B  for  the  spectral  plot. 

Example  3 Two  sinusoids 

Given  = cosn/4.k  + cos.95ir/4.k  for  lags  k = 0,...,8 

Maximum  predictor  order:  M = 4 

Covariance  Sequence  Approximant : 

R^  = (.5000  + j .349910-4)  (1.000  exp  j . 7854)  ^ k I + c.c. 

+ (.5000  - j .349910-4) (1.000  exp  j .7461) ^ + c,c. 

To  4 decimal  accuracy  ir/4  = .7854  and  .95  tt/4  = .7461,  which 
implies  accurate  determination  of  both  closely  spaced  frequencies. 


10  log10|s(z«eJ  ) I (in  dB) 


FREQUENCY  (in  ir  rad) 


Figure  4.2A  Spectral  estimate  for  sine  in  white  noise.  Given  7 known 
covariance  lags.  SNR  - 40  dB. 
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The  corresponding  power  estimates  are  equally  accurate.  See  Figure  4.3 
for  visual  spectral  information. 

Example  4 Two  sinusoids  in  white  noise 

A.  Given  * cosir/4.k  + .01  cos.95ir/4.k  + 100  6(k)  for  lags 

k = 0, ...  ,14  (SNR— 20dB,  -40dB) 

Maximum  predictor  order:  M * 5 

Covariance  Sequence  Approximant: 

-S  I kl 

Rfc  = (.5000  + j -A9121()  J)(1.000  exp  j .7854)'  1 + c.c. 

+ (.49971q“2  - j .491810-5) (1.000  exp  j -7461) ^ + c.c. 

+ <.iooo103  )(.iooo10-10  )'k* 

Note  both  closely  spaced  frequencies,  their  quite  different  powers 
and  the  power  of  the  noise  they  are  submerged  in,  are  all  identified  to 
at  least  3 decimal  places.  See  Figure  4.4A. 

B.  Given  = cosir /4.k  + .001  cos.95ir/4.k  + 1000.  d(k)  for  lags 
k = 0,...,11  (SNR-30dB,  -60dB) 

Maximum  predictor  order:  M * 5 
Covariance  Sequence  Approximant: 

Rj^  - (.4993  - j .575910-3)  (1.000  exp  j .7854)  ^ + c.c. 

+ (.123410-2  + j .576310“3)(.9935  exp  j .7622) I + c.c. 

+ (.10001Q4  )(.103310‘9  exp  j 3.142) ^ 

Single  precision  computation  starts  to  affect  the  least  powerful 
component  here.  Note  that  the  relative  accuracy  is  determined  by  the 
relative  importance  of  each  component  in  terms  of  power.  For  visual 
spectral  information  see  Figure  4.4B,  and  note  that  frequency  .95  ir/4 

-3 

did  not  have  a pole  with  radius  within  1-^q  of  one  to  be  recognized 
as  a pure  sinusoid,  and  the  component  is  therefore  obviously  buried 


in  the  white  noise  component. 


FREQUENCY  (in  it  rad) 


Figure  4.4A 


Spectral  estimate  for  2 sines  in  white  noise.  Given  15 
known  covariance  lags.  SNR  - -20  dB,  -40  dB. 


i 
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Example  5 A sinusoid  in  autoregressive  noise 

A.  Given  Cfc  - cosir /4.k  + c£  for  lags  k - 0 19 

where  is  the  covariance  sequence  for  the  process 
resulting  from  driving 

H(z)  = — r-  4.23 

l-.8z 

with  white  noise  (0,1).  It  is  easily  derived  that 

Cn  = = 2.778 

° 1-.82 

Maximum  predictor  order:  M * 5 
Covariance  Sequence  Approximant: 

\ =■  (.5000  + j .392110-7) (1.000  exp  j .7854) ^ + c.c. 

+ (2.778  ) ( . 8000  )lkl 

The  frequency  estimate,  corresponding  power  and  also  the  first 
order  autoregressive  noise  are  easily  seen  to  be  accurate  to  at  least 
4 decimal  places.  The  visual  spectral  information  is  given  in 
Figure  4.5A.  Note  that  the  maximum  order  was  not  used  here,  since  the 
algorithm  found  accurate  representation  with  order  3. 

B.  As  above  but  now  c”  is  the  covariance  sequence  resulting  from 
white  noise  (0,1)  input  to 

H(z) — 4.24 

(l-.8iz  )(l+.8iz"x) 

The  sequence  was  generated  according  to  the  procedure  outlined 
in  [D].  Covariance  Sequence  Approximant: 

“ (.5000  + j .1688^0  8)(1.000  exp  -j  .7854) ^ + c.c. 

+ (.8469  + j .619510“8)(.8000  exp  -j  1.571) + c.c. 

+ (•143310~7  ) (.7378  exp  -j  3.142) I kl 
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All  estimates  were  found  to  be  accurate  to  at  least  4 decimal 

places.  A spurious  pole  was  Identified  here,  but  the  correction  Is 

achieved  by  finding  an  appropriate  negligible  linear  weight  factor  of 
—8 

order  10  . The  spectral  information  is  depicted  in  Figure  4.5B. 


Example  6 A nonrational  spectral  density 

2 

A.  Given  = exp(-.lk  ) for  lags  k = 0,...,34 

which  represents  the  Gaussian  spectral  density.  Since  it  is 
nonrational,  we  know  there  is  not  any  ARMA(N,M)  that  represents 
it  exactly.  This  example  therefore  serves  as  an  experimental 
test  for  the  robustness  of  the  covariance  sequence  approximant. 
Maximum  predictor  order:  M = 10 
Covariance  Sequence  Approximant: 


= (*3724103  + j .1771103)(.1536  exp  j .7666) ^ + c.c. 

+ (-.31471q3  - j .87521q3)(.1536  exp  j .2503) ^ + c.c. 

+ (.57221q2  - j .42821q2)(.1533  exp  -j  1.348) ^ + c.c. 
Suppose  the  sampled  covariance  sequence  (C^)  is  derived  from  the 


continuous  covariance  function 

2 2 

/ \ ~a  T 
p(x)  = e 

by  sampling  at  = k,  k * 0,±1, 

associated  with  p(t)  is 

. /iT  -(m2/4o2) 
5(0’)  - — e 


4.25 

The  spectral  density  function 

4.26 


With  * k the  sampling  frequency  equals  one  and  the  spectral  density 
of  (4.26)  is  down  107  dB  at  the  half  sampling  frequency  relative  to  the 


zero  frequency.  The  aliasing  effect  in  the  infinite  sura  of  translates 
of  the  continuous  time  specLrum  is  therefore  negligible,  and  we  can  use 


continuous  time  spectral  density  values  in  Table  4.1  to  get  some  idea 
of  the  performance  of  the  covariance  sequence  approximant  procedure. 


w(in  rr  rad) 

10  log10S(w)(in  dB) 

0. 

7.49 

.05 

7.22 

.10 

6.41 

.15 

5.09 

.20 

3.20 

.25 

.79 

.30 

-2.16 

.35 

-5.64 

.40 

-9.66 

.45 

-14.2 

.50 

-19.3 

.60 

-31.1 

.70 

-45.0 

.80 

-61.1 

.90 

-79.3 

1.0 

-99.7 

Table  4.1  Gaussian  spectral  density  values 


The  spectral  information  in  Figure  4.6A  shows  that  the  estimate  is 
fairly  good  up  to  about  oj  = . 5ir  rad.  Sidelobes  then  occur  at  approxi- 
mately -30dB  relative  to  the  peak  spectral  density.  Where  the 
sidelobes  occur  the  spectral  density  estimate  is  negative  for  some 
frequencies. 


Maximum  predictor  oirder:  M * 7 

Note  that  the  major  difference  between  this  example  and  example  6A 

is  the  use  of  about  half  as  many  covariance  lags. 

Covariance  Sequence  Approximant : 

= (4.348  + j .6075)  (.2644  exp  j .9222)  ^ + c.c. 

+ (-.20811q2  - j 2.856) (.2634  exp  j .4498) ^ + c.c. 

+ (-.2002  + j .537110_1)(.2662  exp  j 1.462) ^ + c.c. 

+ (.34321q2  ) (.2631  ) ^ 

+ (.3884i0-4  )(.100010_1°  )'k' 

All  7 poles  were  used,  resulting  in  an  almost  perfect  fit  of  all  the 
given  covariance  sequence  lags.  The  automatically  added  noisepole  got 
an  appropriate  negligible  weight  of  order  10  We  also  note  here  that 
no  negative  spectral  density  estimate  was  encountered  in  this  case. 

Table  4.1  indicates  a very  accurate  spectral  density  for  frequencies  up 
to  u)  = .8n  rad,  after  which  the  cumulative  aliasing  effect  becomes 
noticeable.  See  Figure  4.6B  for  the  spectral  plot. 

If  perfectly  known  covariance  sequences  are  available,  often  only 
a finite  part  is  necessary  to  determine  an  underlying  ARMA(M,M~)  model 
that  implicitly  extends  the  observed  covariance  sequence  for  all  lag 
values.  For  most  examples  considered  so  far,  an  ARMA(M,M  ) process 
or  a limiting  version  of  it,  was  the  underlying  process.  All  results 
were  obtained  with  single  precision  arithmetic.  A nonrational  spectral 
density,  such  as  the  Caussian  one,  can  be  approximated  very  closely  with 
an  estimated  underlying  ARMA(M,M  ).  The  quite  successful  rational 
approximation  points  to  a certain  degree  of  robustness  for  the  covari- 
ance sequence  approximant  and  its  associated  spectral  estimate.  We  also 


note  that  in  none  of  the  examples  for  known  covariance  sequences  was  an 
unstable  ARMA(M,M  ) identified  as  the  approximant . 

4.5  Numerical  Results  for  Estimated  Covariance  Sequences 

It  is  now  logical,  and  very  much  of  practical  significance,  to 
evaluate  the  covariance  sequence  approximant  when  operating  on  a finite 
length  covariance  sequence  estimate  derived  from  a finite  length  realiza- 
tion of  an  observed  process.  In  the  examples  to  follow  realizations  are 
generated  for  several  stochastic  processes,  with  and  without  additive 
noise.  From  these  realizations  a finite  length  covariance  sequence 
estimate  is  determined,  which  will  be  used  to  find  the  covariance 
sequence  approximant  and  corresponding  spectral  information. 

Unless  otherwise  indicated,  a 100  sample  realization  is  generated 
from  the  following  algorithm: 

* sin(ir/4.i)  + S2  cos(frac  x/4.1)  + SN  n^  4.27 

where  S2  represents  the  amplitude  of  a second  frequency  component,  which 
has  a frequency  either  . 95rc/4  or  .6it/4,  and  where  SN  controls  the  power 
of  a (0,1)  white  Gaussian  noise  generator.  The  power  at  the  frequency 

x/4  is  kept  at  1/2,  and  the  power  at  the  second  frequency  is  given  by 

2 2 
1/2  • (S2)  , whereas  the  noisepower  is  given  by  (SN)  . From  the  finite 

length  realization  generated  by  (4.27)  a covariance  sequence  estimate  is 

determined  according  to  the  following  time  averaging  algorithm: 

1 1°«- W 

°k  " 100- 1 kT  *iXi+k!  k * 0 19  4-28 

The  estimator  was  chosen  for  its  property  of  unbiasedness,  The  number  of 
lags  is  a purely  heuristic  decision.  If  in  the  estimator  one  divides 
by  the  total  number  of  data  points,  a biased  covariance  sequence 
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estimator  is  obtained,  and  the  covariance  sequence  approxlmant  for  it 
yields  extra  damping  of  the  pole  locations,  as  is  to  be  expected  from 
the  windowing  effect.  The  maximum  order  in  the  following  examples  was 
arbitrarily  set  at  seven.  In  this  way  we  adhere  to  the  philosophy  that 
one  should  avoid  fitting  the  stochastic  covariance  sequence  estimate 
exactly.  Uith  maximum  order  seven  and  the  use  of  twenty  covariance 
sequence  lag  values  the  fit  is  least  squares. 

The  covariance  sequence  approxlmant  will  form  an  approximation  to 
the  covariance  sequence  estimate,  and  it  therefore  seems  logical  to 
investigate  the  behavior  of  the  covariance  sequence  estimate  derived 
from  a finite  data  record.  Let  us  first  determine  that  behavior  for  a 
finite  observation  of  a pure  sinusoid: 


X±  = A sin(0k+$);  0*0  4.29 

Defining  the  operator  with  the  following  relation 

n i M1! 

- imr  \Vi  4-3° 

and  using  trigonometric  identities,  we  derive  for  the  covariance 
sequence  estimate 

£i  ■ a2cHw 

- A2e”{sin0k+<|>  sin(k+i)0+4>} 

A2  N 

- =Y  ei(cos0i  - cos (2k0+0i+2$) } 

2 2 

■ co80i  - e^{cos2k0+2$»co80i  - sin2k0+2$ *81001} 

* cos0i|"l  - e^{cos2k0+2$}l  + sin0ie^{sin2k0+2$} 


4.31 
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The  covariance  sequence  Is  seen  to  consist  of  two  orthogonal  components . 
The  cosOl  component  converges  to  the  theoretically  expected  covariance 
sequence,  and  the  slndl  component  converges  simultaneously  to  zero.  For 
100  datapolnts  we  can  bound  the  error  terms  by  1 percent  of  the  summed 
value  of  samples  over  half  a period  of  the  double  frequency  26.  As  the 
covariance  sequence  approxlmant  has  the  versatility  to  represent  the 
sln61  component,  we  should  expect  to  see  a good  approximation  to  the 
expression  in  (4.31),  rather  than  a good  approximation  to  the  theoreti- 
cal covariance  sequence  for  the  process  in  (4.29)  based  on  ensemble 
averaging  and  given  by: 

C±  - j A2cos0i  4.32 

Note  that  the  model  underlying  the  Pisarenko  decomposition  cannot 
represent  the  orthogonal  sinOl  component  associated  with  the  covariance 
sequence  estimate  of  (4.31).  The  magnitude  of  the  orthogonal  error 
component  furthermore  can  give  an  indication  of  the  error  term  affecting 
the  theoretically  expected  cos6i  term.  If  the  summation  operator  on 
sin2k6  samples  yields  zero  because  it  operates  over  an  integer  number 
of  periods,  then  the  same  summation  operator  on  cos2k6  samples  also 
yields  zero.  In  the  latter  case  the  magnitude  of  the  cosSi  component 
gives  an  indication  of  the  power  in  the  sinusoid. 
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A single  sinusoid 


S2  - SN  - 0 


Covariance  Sequence  Approxlmant: 


Rj^  = (.897210“5  - j .439210"3)  (1.034  exp  j .7851)  ^ + c.c. 

+ (.2500  + j .20611q"2) (1.000  exp  -j  .7854) ^ + c.c. 

+ (-.59951Q_9  + j . 30251Q'9)(. 7555  exp  j 2.6Q8)lkl  + c.c. 


+ (.4079, 


+ (-.3161, 


M.3151  exp  j 3.142)1 


)(.1000, 


As  indicated  in  (4.31)  an  orthogonal  component  could  be  expected 
with  magnitude  less  than  1 percent  of  the  true  power,  which  is  seen  to 
be  valid.  The  frequency  estimate  is  seen  to  be  accurate  to  4 decimal 
places.  Several  spurious  poles  were  identified  but  their  respective 
linear  weights  are  very  small.  Note  the  occurrence  of  an  unstable  pole 
close  to  the  dominant  frequency,  which  is  then  discounted  by  a very 
small  linear  weight  coefficient.  The  spectral  information  depicted  in 
Figure  4.7  was  computed  according  to  the  algorithmic  expression  for  the 
parametric  spectrum  estimate  in  (4.2).  A spectrum  does  not  exist  in 
this  case  because  {R^}  is  explosive  by  virtue  of  the  poles  with  radius 
1.034.  For  an  unstable  pole  pair  however  that  computation  yields  the 
negative  of  the  spectral  density  for  the  stable  pole  pair  associated 
with  the  reflection  of  the  unstable  poles.  For  this  particular  covari- 
ance sequence  approxlmant  the  unstable  pole  pair  forms  the  only  major 
contribution  to  a "continuous  spectral  density"  which  results  in  the 
deceptively  regular  appearance  of  Figure  4.7. 

The  next  question  to  be  answered  is  what  happens  when  more 


sinusoidal  components  are  present?  To  provide  at  least  a partial 


answer  we  examine  the  following  realization: 

Xi  = AjSinOji-h^)  + A2sin(021442)  4.33 

Repeated  use  of  trigonometric  Identities  yields  the  following  covariance 
sequence  estimate: 

1 M1! 

Ci  " N=flf  k£1  (Vin(eik+*l)+A28in(02k^2))(Al8in(0l(fcfi)'H,,l) 

+ A2sin(62(k+i)+*2)) 

= A^e^{sin(01k+<(>1)sin(ei(k+i)-l-<Ji1)}  + A2e^{sin(09k+42)sin(02(k+i)+<|>2)  } 

+ A^A2E^{sin(0^k+4^)sin(02(k+i)+4>2)  + sin(02k+4>2)sin(01(k+i)+$1)} 

= | Ajcos01i[l-eJ{cos(201k+2*1)}J  +|  A^sinO^e^ 81^20^+2^)} 

+ | A2cos02i[l-e“{cos(202k+2<|.2)}]  +|  A2sin02ieJ{sin(202k+2*2) } 

+ \ A^cosOj^i  ^{oos(02-01)k-H|i2-iJ»1}  - e^{cos(01+O2)k+4>1+«|i2}J 
+ i A1A2sin01l  [e^{sin(02-01)k+^2-<|.1}  + e“{sin(01+02)k+*1+*2}] 

+ 7 Aia2cos021  |e^{cos(0j-02)k+<f>j-$2}  - e^{cos(01+02)lcH1+42>J 

+ 1 A1A2sine2i  [ei^sin^ei_02^k+<^i"<*,2^  + ei^sin^0l+02^k+^l+^2^] 

* cosO^y  A^(l-e^{cos201k+2<|i1})  A^e^cos^-O^fc^-^ 

- cos(01+02)k+<t>1+(}i2}J 

+ slnO^  Ajc^{sln201k+2*1}+|  A^cjlsin^-O^fc*^-^ 

+ sin(O1+02)k-H|.1-h|>2}] 

+ cos02i^  A2'(l-e^{cos202k+2t2}  + y A1A2eitcoa(0i“02)k+*i_*2 

- cos(01+02)k+4>1+^2)J 


+ sin02i[i  A^eJ{sin202k+2*2>+i  AiVi^^VV^rt 

+ 8in(01+02)k+41+42>]  4*34 

Note  again  that  a decomposition  into  orthogonal  components  takes  place 
and  furthermore  that  all  error  terms  involving  • } converge  to  zero, 
if  0^02^0.  Thus  the  covariance  sequence  estimate  based  on  a single 
realization  is  an  asymptotically  unbiased  estimator  for  the  covariance 


sequence. 


For  the  examples  to  follow  a process  realization  was  generated 
with  the  algorithm  of  (4.27)  so  that  in  (4.33): 


\ - 1 , a2  * s2 

0^^  = tt/4,  ©2  - frac  if/4 
4>1  ■ 0 . *2  “ 
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Example  8 Two  sinusoids 


A.  S2  = 1 


frac  « .95 


SN  = 0 


Covariance  Sequence  Approximant : 

\ = (.3611  - j .1508) (.9989  exp  j .7819) ^ + c.c. 

+ (.3581  - j .1535) (.9988  exp  -j  .7497)'k'  + c.c. 
+ (.4111.  "5  w ^kl 


) (.8340 


+ (.132110-5  - j .39821q"6)(.5592  exp  j 2.443) M + c.c. 


+ (-.5151, 


H.1000, 


Both  frequency  estimates  are  correct  to  within  1/2  percent.  One 


might  expect  that  this  leads  to  a very  good  overall  estimate,  and  indeed 


all  twenty  lags  are  matched  to  at  least  4 decimal  places.  Spurious 


poles  have  indeed  an  appropriate  low  linear  weighting  coefficient.  The 


power  estimates  however  are  substantially  too  high.  From  (4.34)  we  see 


that  a double  frequency  and  sum  frequency  component  enters  as  error 


term,  and  by  the  same  argument  as  before  they  are  of  the  order  of  1 per- 


cent. The  more  important  error  term  in  this  case  could  be  the  differ- 


ence frequency  term.  This  difference  frequency  is  ir/4  - .95ir/4  = n/80. 


With  a 100  point  data  sequence  this  turns  out  to  be  an  almost  worst  case 


situation  where  the  error  terms  are  very  large.  The  covariance  sequence 


approximant  above  merely  satisfies  the  covariance  sequence  estimate. 


The  spectral  information  is  depicted  in  Figure  4.8A.  Note  that  there  is 


no  spectral  resolution  readily  available  in  the  plot,  whereas  the 


covariance  sequence  approximant  clearly  indicates  two  strong  sinusoidal 


components . 
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FREQUENCY  (in  tt  rad) 


Figure  4.8A  Spectral  estimate  for  2 sines.  Estimated  20  covariance  lags 
from  100  samples. 


i 
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B.  Based  on  the  frequency  estimates  from  example  8A  and  the  error 
terms  indicated  in  (4.34),  we  can  decide  to  minimize  these  error  terms 
by  insuring  that  summation  takes  place  over  an  integer  number  of 
periods  of  the  difference  frequency.  Let  us  therefore  change  the 
length  of  the  data  sequence  to  N - 160,  with  all  other  parameters 
unchanged. 

Covariance  Sequence  Approximant: 

Rj^  = (-.103210-8  - j .470010'9)(.9599  exp  j 2.236)  ^ + c.c. 

+ (.2500  - j .170810-2) (1.000  exp  -j  .7462) + c.c. 

+ (.2500  - j .154510“2) (1.000  exp  -j  .7854) ^ + c.c. 

+ (-.377010-8  ) (. 7687  exp  j 3.142) 

+ (-.100510"7  )(.100010"10  )lkl 

Both  frequencies  are  very  accurate,  and  also  their  respective 
power  estimates  are  accurate  to  4 decimal  places.  Spurious  poles  were 
assigned  appropriate  low  weights,  and  the  orthogonal  components  are 
of  the  order  of  = .00625  percent  as  is  to  be  expected.  See 
Figure  4.8B.  This  example  shows  therefore  that  knowing  the  component 
frequencies,  the  observation  interval  can  be  chosen  to  eliminate  the 
error  components  due  to  difference  frequencies.  Such  a procedure  is 
of  course  not  possible  if  the  availabK  data  record  is  shorter  than  one 
period  of  the  difference  frequency  of  interest. 

Let  us  now  assume  the  ultimate  practical  case  where  our  observation 
interval  is  finite  and  the  observations  are  affected  by  additive  noise. 

For  white  noise  we  then  make  an  estimate  of  the  covariance  sequence  for 
lags  different  from  zero: 


CD 

CM 

I 


0.00 


"0V2O 


0.40  0.60  0.80 

FREQUENCY  (in  n rad) 


1 . 00 


Figure  4.8B  Spectral  estimate  for  2 sines.  Estimated  20  covariance  lags 
from  160  samples. 


is  white  noise,  so  that 


\ " W+i 


As  a result  the  white  noise  covariance  estimate  at  lag  i j*  0 has  a 


standard  deviation 


Since  the  white  noise  process  is  deemed  independent  from  the  signal 
process,  this  result  is  additive  to  results  derived  previously  for  co- 
variance  sequence  estimates  from  finite  record  length  noiseless  data. 
This  elementary  calculation  simply  illustrates  that  large  variance  white 
noise  processes  observed  over  short  data  records  yield  large  variance 
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Example  9 A single  sinusoid  In  white  noise 

A.  S2  - 0 

SN  - .03  (SNR  - 27dB) 

Covariance  Sequence  Approxlmant: 

Rfc  - (.2481  - j ,247210-2) (1.000  expj  .7850)lkl  + c.c. 

+ (.710310“3  ) (.7826  )lkl 

+ (.712210‘3  + j .209810_4)(.6216  exp  -j  2.668) ^ + c.c. 

+ (.413210-2  - j .766310‘3)(.3997  exp  -j  .9054) ^ + c.c. 

+ (-.779510-2  )(.100010-10  )lkl 

The  sinusoidal  pole  and  its  corresponding  power  are  identified 
quite  nicely.  Spurious  poles  are  identified  and  seem  to  be  approximating 
the  noise  process.  See  Figure  4.9A  for  visual  spectral  information. 

B.  S2  = 0 

SN  = .3  (SNR  = 7dB) 

Covariance  Sequence  Approxlmant: 

= (.2364  - j .247010-2)  (1.000  exp  j .7827)  I kl  + c.c. 

+ (.202610_1  ) (.9295  )lkl 

+ (-.128510_1  - j .82711q"2)(.7551  exp  -j  2.005) • + c.c. 

+ (.1335  + j .1051) (.3341  exp  -j  2.206) M + c.c. 

+ (-.1614  )(.100010~10  )'kl 

NoLc  the  frequency  estimate  is  very  accurate  and  the  corresponding 
power  estimate  is  reasonable.  Spurious  poles  are  such  that  they  repre- 
sent an  average  noise  level  of  -lOdB  as  well  as  can  be  expected.  See 
Figure  4.9B. 


■' — - 


10  l°g10 I S(z=e^  9) | (in  dB) 
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Figure  4.9A  Spectral  estimate  for  sine  in  white  noise.  SNR  * 27  dB. 
Estimated  20  covariance  lags  from  100  samples. 
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C.  s2  = 0 

SN  - .3  (SNR  = -13dB) 

Covariance  Sequence  Approxlmant : 

Rfc  = (.5260  - j 2.549) (. 7368  exp  j 1.255) ^ + c.c. 

+ (3.530  - j .4836) (. 6648  exp  -j  1.661) + c.c. 

+ (1.093  ) (. 9881  )!kl 

+ (1.278  + j . 7616) (. 7755  exp  j 2.250) ^kl  + c.c. 

-(2.138  )(.100010'10)lkl 

As  indicated  before  the  stochastic  nature  of  the  noise  covariance 
estimate  might  take  importance  over  any  signal  features  when  the  noise 
variance  is  comparatively  large.  The  identified  poles  therefore 
represent  the  noise  process  and  as  seen  in  Figure  4.9C  the  spectral 
density  estimate  varies  around  the  lOdB  level. 

Example  10  Two  sinusoids  in  white  noise 
A.  S2  = 1 

SN  = .03  (SNR  = 27  dB,  27dB) 
frac  = .95 

Covariance  Sequence  Approxlmant: 

Rk  - (.8815  + j .796310-4)(.9973  exp  j ,7655) lk|  + c.c. 

+ (.54771q~3  ) (.5317  exp  j 3.142)lkl 

+ (-.1608  + j .223310"2)(.9605  exp  j ,7630) lkl  + c.c. 

+ (-.569410-5  + j .28181q'3)(.6813  exp  -j  2.226) lk>  + c.c. 
+ (.  7627jQ-3  )(.10001()-10  )lkl 


4 


Figure  4.9C  Spectral  estimate  for  sine  In  white  noise.  SNR  - -13  dB. 
Estimated  20  covariance  lags  from  100  samples. 
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No  resolution  of  the  frequency  components  is  obtained.  Instead  the 
frequencies  are  virtually  averaged,  which  may  then  cause  a large  devi- 
ation in  the  respective  power  estimates  due  to  error  terms  dependent  on 
the  difference  frequency.  In  this  case  one  of  the  frequencies  shows  a 
considerable  negative  power  estimate.  The  spectral  information  is 
depicted  in  Figure  4.10A.  Note  from  (4.39)  that  the  stochastic  element 

— A 

of  the  covariance  sequence  elements  is  of  the  order  10  * which  is 

substantial  compared  to  the  difference  between 

C2  - y cosir/4.k  + ^ cos.95ir/4.k 

and 

C3  * co8.975ir/4.k 
over  the  lag  values  k = 0,...,19. 

B.  S2  - 1 

SN  - .03  (SNR  3 27dB,  27dB) 
frac  ■ .6 

Covariance  Sequence  Approximant: 

Rk  - (.2486  + j .351310-2) (1.000  exp  j .4717) M + c.c, 

+ (.100310-3  - j .223410"3)(.8381  exp  j 2.121) lkl  + c.c. 

+ (.2492  - j .335410*2) (1.000  exp  -j  .7848) + c.c. 

+ (.263310“3  ) (.7696  exp  -j  3.142) ^ k ^ 

+ (.181510“3  )(.100010‘10  )lkl 

Both  frequency  estimates  are  within  .1  percent  and  their  power 
estimates  are  within  1 percent.  Note  that  the  difference  frequency  is 
ir/10  so  that  with  100  datapolnts  the  error  terms  due  to  finite  record 
estimates  are  vanishingly  small.  See  Figure  4.10B  for  spectral 
Information  plot. 


*iols< 


k 

B 


109 

c.  s2  - 1 

SN  - .3  (SNR  • 7dB,  7dB) 
frac  - .6 

Covariance  Sequence  Approximant : 

= (.2589  - j .167610'1)(.9940  exp  j .4722) ^ + c.c. 

+ (.491610-3  ) (. 7588  exp  j 3.142) I kl 

+ (.2442  + j .258410_1)(.9997  exp  j . 7749) + c.c. 

+ (.262510‘1  - j .526910"3)(.5969  exp  -j  2.051)^  + c.c. 

+ (-.125410_1  )(.100010~10  )lkl 

Power  and  frequency  estimates  are  still  within  3 percent  and  the 
extra  poles  seem  to  represent  a nolselevel  of  -lOdB  rather  well  In 
Figure  4. 10C. 

D.  S2  = 1 

SN  = 3 (SNR  * -13dB,  -13dB) 
frac  = .6 

Covariance  Sequence  Approximant: 

Rk  = (-.2295  - j .4725) (.9087  exp  j 1.311)'k'  + c.c. 

+ (-.859310-2  - j .4811) (.8643  exp  -j  2.122) M + c.c. 

+ (.8077  ) (1.011  )lkl 

+ (.5078101  - j .5756101)(.3821  exp  j 2.245) ^ + c.c. 

+ (-.6262  )(.100010‘10  )lkl 

The  sinusoidal  frequency  components  are  burled  in  the  noise  and  can- 


not be  recovered  even  in  an  averaged  frequency  component.  The  stochastic 
element  of  the  covariance  sequence  estimate  according  to  (4.39)  is  now 


Ill 


of  the  order  1 which  is  the  maximum  value  of  the  covariance  sequence  for 
the  averaged  frequency.  Note  in  the  plot  of  Figure  4.10D  the  attempt 
to  represent  a noise  level  of  lOdB.  Note  again  the  unstable  pole,  so 
that  Figure  4.10D  really  reflects  the  spectral  density  as  if  the  linear 
weight  for  the  unstable  pole  were  negative  and  associated  with  the 
reflection  of  the  unstable  pole. 

Finite  record  based  covariance  sequence  estimates  can  drastically 
affect  the  power  estimates  corresponding  to  sinusoidal  frequencies.  A 
one-shot  procedure  is  therefore  not  likely  to  produce  simultaneously 
good  frequency  and  power  estimates  when  analyzing  short  data  records. 
Based  upon  the  frequency  estimates  that  result  from  the  first  step  in 
the  determination  of  the  covariance  sequence  approximant,  the  most 
favorable  length  of  the  analysis  record  may  be  determined.  Calculating 
the  complex  linear  weight  coefficients  for  the  covariance  sequence 
approximant  of  this  thesis  from  this  adjusted  analysis  record  will  mini- 
mize the  error  components  due  to  finite  length  observation  records,  as 
in  (4.34).  Uhen  observations  are  in  addition  noisy,  the  variance  of 
the  noise  covariance  estimate  may  obscure  the  signal  covariance  estimate. 
The  result  is  a loss  of  resolution.  Initially  a single  frequency  com- 
ponent is  identified  at  averaged  frequency  and  power  of  two  separate 
components.  For  even  noisier  data  an  approximation  of  the  noise  process 
takes  place. 


) . 60  0.80 

FREQUENCY  (in  tt  rad) 


Spectral  estimate  for  2 sines  in  white  noise.  SNR  = - 
-13  dB.  Estimated  20  covariance  lags  from  100  samples 
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5 COVARIANCE  SEQUENCE  APPROXIMATION  FOR  RECURSIVE 


DIGITAL  FILTER  DESIGN 

In  the  previous  chapter  we  have  seen  that  the  covariance  sequence 
approxlmant  leads  to  excellent  spectral  analysis  results  when  it  is 
applied  to  known  covariance  sequences.  As  the  model  underlying  the 
approxlmant  is  an  ARMA(M,M~)  system,  it  seems  natural  to  use  the 
covariance  sequence  approxlmant  to  design  ARMA  recursive  digital 
filters,  based  exclusively  on  covariance  sequence  information.  In  the 
first  section  of  this  chapter  we  restate  a few  properties  derived 
earlier  and  indicate  how  the  possibly  negative  spectrum  estimate 
associated  with  the  covariance  sequence  approxlmant,  can  be  avoided 
in  a slightly  modified  covariance  sequence  approximation  procedure 
for  the  design  of  recursive  digital  filters.  The  resulting  filter 
can  be  guaranteed  to  be  stable,  causal,  and  have  real  coefficients. 

In  Section  5.1  we  describe  the  freedom  available  in  the  choice  of 
phase  characteristics  that  can  be  associated  with  the  resulting  magni- 
tude squared  approximation.  This  relative  freedom  for  the  phase 
characteristic  is  associated  with  the  different  possibilities  for  the 
numerator  polynomial  of  the  filL^r  and  several  ways  are  indicated  to 
arrive  at  a possible  set  of  coefficients.  Section  5.3  contains  a 
compilation  of  numerical  considerations  associated  with  the  design 
procedure  of  this  thesis  and  the  final  section  contains  results  for 
several  representative  low  pass  recursive  digital  filter  designs. 

5.1  Recursive  Digital  Filter  Design 

Suppose  we  wish  to  design  a recursive  digital  filter  H(z)  so  that 
its  associated  spectrum  forms  a good  approximation  to  a desired 
spectrum.  Many  design  procedures  have  been  proposed,  starting  from 
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Che  frequency  domain  [D] , [KC],  [DU],  or  from  Che  Cime  domain  [BS], 

[EF] . In  general  Che  associaCed  lease  squares  problem  is  highly 
nonlinear,  which  forces  upon  us  a Cime  consuming  iceracive  mode  of 
soluCion.  Simple  procedures  based  on  impulse  response  approximaCion 
have  generally  noC  succeeded  Co  guaranCee  sCabiliCy  of  Che  resulcing 
filCer  [BP],  [SH].  To  arrive  aC  a more  amenable  problem  sCaCemenC, 
modified  lease  squares  procedures  have  been  proposed  [SM],  [MI]  and 
chis  has  recenCly  led  Co  an  efficienC  design  procedure  for  sCable  ARMA 
(N,M)  digiCal  filcers  from  consisCenC  impulse  and  covariance  sequence 
data  [MR] . Design  resulCs  wich  Che  laCCer  procedure  have  been  reporced 
[SL]  and  clearly  show  Che  influence  of  Che  modified  lease  squares 
error  criCeriou. 

We  have  noced  before  ehac  an  ARMA(M,M  ) model  was  Che  appropriace 
structure  for  modeling  when  it  is  desired  to  preserve  second  order 
statistical  properties  in  transforming  continuous  time  rational  pro- 
cesses to  discrete  time  processes  [PS].  An  M*"^1  order  continuous  time 
process  leads  to  an  ARMA(M,M  ) discrete  time  process,  such  that 

C(k)  = Cc(x-kT)  Vk  5.1 

A covariance  sequence  so  derived  exactly  satisfies  an  Mth  order  linear 
homogeneous  differencial  equation  with  real  coefficients.  The  covari- 
ance sequence  approximant  yields  zero  covariance  sequence  prediction 
error  and  Che  resulCs  are  exact.  The  covariance  sequence  approximaCion 
procedure  Chen  merely  serves  as  an  alternative  procedure  for  Cl 
(covariance  invariant)  digital  filter  design  [PS],[KPSS].  The  more 
difficult  task  arises  when  there  is  no  ARMA(M,M  ) or  oven  AHMA(N.M) 
system  that  matches  {Ck}QN  1 unless  N=2M.  The  Latter  condition  would 
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constitute  a Padd-like  approach  to  covariance  sequence  approximation 
and  in  general  would  lead  to  a much  higher  order  for  the  approximating 
system  than  desirable.  Also  In  this  case  nothing  would  be  known  about 
the  match  of  {C^}  for  lags  larger  than  N. 

We  will  now  indicate  how  a digital  filter  with  real  coefficients 
can  be  designed,  using  the  covariance  sequence  approximant  as  derived. 
Let  Sd(0),  -it<6<it,  denote  an  ideal  spectrum,  with  corresponding  co- 
variance  sequence  (c£)  e i2  , to  be  approximated  with  an  ARMA(M,M~) 
digital  filter.  The  least  squares  approximation  problem  is  to  find 
the  parameters  {A^,p^}^  for  the  covariance  sequence  approximant  e 


5.2 


5.3 


is  minimized.  Here  Sd(0)  is  the  spectrum  estimate  corresponding  to 
the  ARMA(M,M  ) covariance  sequence  approximant.  Note  the  fundamental 
discrepancy  between  the  %2  sequence  to  be  approximated  and  the 
approximating  sequence.  The  problem  stated  above  is  exactly  the 

general  nonlinear  problem  of  (5.2)  that  we  seek  to  avoid.  The  modified 


version  of  this  problem  has  been  outlined  in  Chapter  3.  The  main 
difficulty  encountered  in  the  direct  application  of  the  modified  least 

A 

squares  problem  is  that  8^(0)  ®sy  turn  out  to  be  negative  for  some 
frequencies,  implying  that  the  sequence  {R^}  is  not  a covariance 

A 

sequence  and,  furthermore,  that  a factorization  of  sj(z)  to  obtain 
the  corresponding  filter  will  give  a complex  filter.  The  negativity 
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of  the  spectrum  estimate  is  somewhat  tolerable  In  the  spectral  analysis 
case,  as  the  problem  Invariably  turns  up  in  the  approximation  of  the 
low  energy  parts  of  the  spectrum.  In  digital  filter  design  it  is 
imperative  that  a filter  with  real  coefficients  result,  and  when  de- 
signing ideal  filters  with  very  low  energy  stopbands,  zeros  on  the  unit 
circle  cannot  be  avoided.  We  note  here  that  in  Cl  design  an  ARMA(M,M  ) 
system  with  real  coefficients  is  known  to  exist  and  will  be  identified. 


but  in  the  general  design  problem  there  is  no  guarantee  that  zeros  of 

S.(z)  on  the  unit  circle  will  occur  in  multiples  of  two.  To  avoid  the 
a 

negativity  problem  we  instead  approximate  the  covariance  sequence 

— 

associated  with  S,  (6),  which  results  in  a ARMA(M,M  ) covariance 
d 

* k 

sequence  approximant  with  associated  spectral  density  estimate  (6). 
We  now  call 

sd(e)  - (Sj^e))2  5.4 

the  approximation  to  S^(6).  Substitution  of  (4.3a)  yields: 


Sd(z) 


A(z)  ACz”1)  A(z  1)  A(z) 

Jsisili 

|a(z) |4 


5.5 


Note  that  S.(z)  satisfies  all  requirements  for  a proper  spectrum, 
a 

The  properties  of  reality  and  symmetry  are  obviously  preserved.  It  has 

the  same  poles  and  zeros  as  S,  (z),  but  all  occur  in  multiples  of  two. 

a 

A 

The  important  aspect  is  that  zeros  of  S,(z)  on  the  unit  circle  now 

a 

occur  in  multiples  of  two  so  that  the  coefficients  of  the  associated 

_ j 

ARMA(2M,2M  ) filter  are  real. 

We  have  seen  that  the  covariance  sequence  approximation  algorithm 
nnn  be  implemented  easily  to  yield  the  stable  denominator  polynomial 
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A(z)  of  maximum  order  M for  a given  covariance  sequence.  Also  the 
coefficients  of  A(z)  are  always  real.  We  note  that  a change  In  length 
of  the  covariance  sequence  to  be  approximated  may  affect  the  order  as 
well  as  the  performance  of  the  resulting  covariance  sequence  prediction 
filter.  To  preserve  stability  the  obvious  choice  for  the  denominator 
polynomial  of  the  recursive  digital  filter  with  spectrum  S^(z)  will 
be  A(z)  A(z).  The  numerator  polynomial  can  be  arranged  in  many  differ- 
ent ways  as  shown  In  the  following  section. 

5.2  Numerator  Polynomial  Coefficients  and  Phase  Characteristics 

We  are  Interested  In  the  zeros  of  S^(z)  in  (5.5).  As  they  occur 

i* 

in  multiples  of  two  It  suffices  to  study  the  zeros  of  Sj  (z).  After 

d 

finding  the  zeros  of  (z)  or  equivalently  the  roots  of  N(z),  the 

A 

complex  roots  of  (z)  with  radius  different  from  one  can  be 
arranged  In  sets  as  follows: 


Four  of  these  zeros  can  be  chosen  for  the  ABMA(2M,2M~)  system,  provided 
they  occur  in  complex  conjugate  pairs,  so  as  to  Insure  the  recursive 
digital  filter  has  real  coefficients.  Seal  zeros  away  from  the  unit 
circle  occur  In  following  sets: 


, -1  -1-.  - 

{po’  po  » V po  } 5 

Two  of  these  roots  are  to  be  chosen,  with  no  restrictions.  If  real 

roots  occur  on  the  unit  circle,  one  of  each  multiple  of  two  Is  to  be 

taken.  Complex  zeros  on  the  unit  circle  occur  In  foursomes 


, 10  -18  10  -10, 


5.6c 


As  they  again  have  to  be  taken  in  complex  conjugate  pairs,  there  is  only 
one  way  to  do  so.  As  a result  of.  the  above,  the  numerator  polynomial 
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of  the  digital  filter  has  real  coefficients.  Recall  that  with  a 
denominator  polynomial  A(z)  A(z)  then  all  digital  filter  coefficients 
will  be  real. 

Note  that  for  the  zeros  away  from  the  unit  circle  several 
combinations  of  zero  pairs  will  be  possible  and  this  allows  considerable 
freedom  in  the  phase  characteristic  to  be  associated  with  the  ARMA 
(2M,2M  ) filter.  It  is  important  to  note  that  this  design  freedom  is 
implicit,  as  it  was  the  magnitude  squared  or  spectral  performance  only 
that  prompted  this  design  procedure.  We  may  therefore  achieve  approxi- 
mately linear  phase  response  over  a frequency  region  of  importance  by 
cleverly  choosing  zeros.  Alternatively  one  may  achieve  partial  phase 
equalization.  Furthermore,  if  the  zeros  can  be  chosen  to  lie  inside 
the  unit  circle,  an  invertible  digital  filter  will  have  been  designed. 


To  exercise  the  available  freedom  of  phase  characteristic  an  additional 
polynomial  rootfinding  will  have  to  be  executed,  even  if  only  on  a 
polynomial  of  order  M~. 

A straightforward  procedure  to  arrive  at  a possible  numerator 
polynomial  is  to  recognize  that  the  following  system  is  causal,  has 
real  coefficients  and  exhibits  spectrum  S^(z)  when  ..riven  by  white 
noise: 


H(z) 


A(z)  A(z) 


k > M~ 


5.7 


where  again,  as  in  (4.3b): 


M M 

N(z)  « E A (1-p.Ml+p.)  n (z-p.Hz"1-?.  ) 

i-1  1 1 k k 


; -j 

- z _ n z ; n 

J— M J 3 


are  real,  n^  - n_^ 


5.8 
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The  coefficients  of  N(z)  can  be  determined  In  straightforward  fashion, 
without  polynomial  root finding  or  spectral  factorization  techniques. 

One  may  even  use  to  advantage  the  fact  that  N(z)  is  a mirror  image 
polynomial,  which  cuts  the  required  computational  effort  roughly  in  half. 

As  the  latter  procedure  still  may  seem  to  be  an  exhaustive  approach, 
s more  appealing  method  will  be  outlined  next.  Recall 


«A.)  - — ! *-&■ 


A(z)  A(z  1) 


M 

- Z A. 


(1-Pi)(l+P1) 


i-1  1 (z-Pi)(zrl-Pl) 


5.9 


Por  any  z - e ^ it  is  therefore  possible  to  evaluate  N(z  - e^0) 
according  to 


N(z  - eJ6)  - A(z)  A(z"1)  §d*(z) I j0 

| z ■ eJ 

- |A(z)|2  S.\z)  ,,  5-10 

z ■ eJ 

Note  that  all  elements  of  the  right  hand  side  of  (5.10)  are  readily 
available  as  a^,  p^  and  A^  have  all  been  determined.  Furthermore 
N(z)  is  real  and  symmetric  on  as  is  easily  seen. 

With  (5.8)  and  (5.10)  we  may  therefore  write 

nt  - IDFT  1 1 A(z)  |2  S^U)!  . #J2*k/2M|  5.11 


The  latter  approach  can  also  be  used  to  evaluate  the  coefficients 
of  N(z)  N(z  ) in  (5.5).  As  N(z)  represents  a pure  moving  average 
system  with  real  coefficients,  the  resulting  sequence  would  represent 
the  covariance  sequence  associated  with  N(z).  Knowing  however  a 
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possible  set  of  coefficients  for  the  numerator  polynomial,  say  from 
(S . 11) , the  corresponding  covariance  sequence  can  be  readily  calculated . 
M-l 

1 nt  n.+v  • k " -2 (M-l) 2 (M-l) 

A- (M-l)  5.12 

0 otherwise 

Once  the  moving  average  covariance  sequence  is  known,  in  addition  to  the 
order  of  the  generating  system,  algorithms  are  available  to  determine 
the  corresponding  minimum  phase  numerator  polynomial  for  the  recursive 
digital  system  by  factorization  techniques  [Wl]. 

3 . 3 Numerical  Considerations 

In  the  following  section  we  will  design  several  low  pass  filters. 

\ * 

The  starting  point  is  naturally  the  ideal  low  pass  spectral  density, 
but  a modification  will  be  made  in  that  the  stopband  level  is  not  speci- 
fied to  be  zero.  Instead  the  stopband  level  is  thought  to  represent 
a hypothetical  tolerable  noiselevel.  We  note  here  that  minimization 
of  the  covariance  sequence  prediction  error  represents  a spectral 
whitening  approach.  The  numerical  conditioning  of  the  system  of  linear 
equations  to  be  solved  for  minimization  of  the  covariance  sequence  pre- 
diction error  is  determined  to  a large  extent  by  the  dynamic  range  of 
the  desired  spectrum  [MG].  Specification  of  a tolerable  noiselevel  for 
the  stopband  performance  of  the  modified  ideal  filter  will  therefore 
stabilize  the  numerical  behavior  oT  the  algorithm. 

Use  of  the  minimal  length  covariance  sequence  to  find  an  approximate 
ARMA(M,M  ) system  leads  to  a Pade-like  approximant  to  the  covariance 
sequence.  In  the  application  of  complex  exponential  approximation  to 
one-sided  impulse  response  sequences  similar  design  procedures  are 
known  [BP],  (SliJ,  [PKj,  [HlJ  and  several  possible  problems  have  been 
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reported  [VB] . Complex  roots  not  occurring  In  complex  conjugate  pairs 
Indicate  that  the  rootfinding  routine  has  not  converged  yet,  as  the 
polynomial  itself  is  known  to  have  real  coefficients.  The  obvious 
solution  to  such  a problem  is  to  allow  more  iterations  in  the  polynomial 
routine.  In  all  covariance  sequence  approximations  performed  in  this 
thesis  the  polynomial  roots  always  were  either  real  or  occurred  in 
complex  conjugate  pairs  within  the  specified  numerical  accuracy. 

Partial  realization  procedures  and  Pad£  approximation  are  known  to 
be  very  sensitive  to  noise  [VB],  [W2].  Host  algorithms  solving  the 
partial  realization  problem  are  for  that  reason  numerically  unstable,  but 
numerically  stable  algorithms  have  been  proposed  requiring  a slightly 
increased  number  of  operations  [DJ].  In  this  respect  we  note  that  in  the 
covariance  sequence  approximation  procedure  of  this  thesis,  equations 
are  always  solved  in  a least  squares  sense.  The  order  of  the  approx- 
imating system  is  low  compared  to  the  length  of  the  available  covar- 
iance sequence.  This  approach  has  a noise  reducing  and  numerically 
stabilizing  effect.  In  the  proposed  digital  filter  design  procedure 
this  sensitivity  problem,  which  would  lead  to  a spectral  factorization 
with  complex  coefficients,  is  avoided  entirely  by  design.  Noise 
sensitivity  does  not  play  the  role  it  does  in  partial  realization 
procedures  because  of  the  fact  that  a covariance  sequence  approximant 
for  the  intermediate  covariance  sequence  does  not  need  to  have  the 
nonnegative  definiteness  property.  The  final  result  for  the  spectral 
approximation  is  still  guaranteed  to  be  nonnegative. 


Ideal  filters  inevitably  have  discontinuous  transitions  in  the 
frequency  domain,  whereas  their  corresponding  covariance  sequences  will 
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be  relatively  smooth  and  damped.  The  approximation  of  such  a two-sided 
covariance  sequence  with  a weighted  sum  of  complex  exponentials  naLurally 
avoids  singular  cases  of  approximating  time  domain  discontinuities 
which  reportedly  yields  serious  problems  [MCD] . The  modified  low  pass 
filters  have  a square  summable  covariance  sequence  whereas  our 
covariance  sequence  approxlmant  is  absolutely  summable,  which  points  to 
a fundamental  mismatch  between  approximated  and  approximating  sequence. 

In  the  approximation  of  known  covariance  sequences  for  both  the 
spectral  estimation  simulations  of  Chapter  4 and  the  digital  filter 
designs  of  the  next  section,  the  denominator  polynomials  are  always 
minimum  phase  corresponding  to  stable  digital  filters.  It  is  even  more 
noteworthy  that  these  results  have  been  obtained  without  recursively 
checking  for  stability  in  order  to  end  up  with  a reduced  order  but 
guaranteed  stable  filter. 

5.4  Numerical  Results  for  Recursive  Digital  Filter  Design 

Let  us  assume  a normalized  sampling  frequency  of  1 and  a normalized 
cutoff  frequency  of  .1  so  that  the  spectral  density  specification  in 
Figure  5.1  is  relevant.  A similar  starting  point  arises  in  the  design 
of  digital  filters  from  a statistical  point  of  view  [FS],  [SL],  The 
known  covariance  sequence  associated  with  the  above  desired  spectral 
density  is  derived  by  sampling  the  Fourier  transform  associated  with 
Sd(f>! 

Cfc  • 2W  slnc(2nWk)  + All  cos{(l+2W)  ^4  (1-2W)  sino{(l-2W)  5.13 

From  Fourier  transform  properties  we  may  intuitively  say  that  the  use 
of  long  sequences  {c^},  will  be  reflected  in  emphasis  on  low  frequency 
approx  I mu  I I on . The  iimc  ol  short  sequencer!  will  cor  respond  to  relatively 
more  emphasis  on  the  high  Irequeney  approximation  error.  When  minimal 


length  sequences  are  used  a Pade  approximation  to  the  covariance 
sequence  results,  with  emphasis  on  the  high  frequency  approximation 
error  and  usually  characteristic  peaking  at  the  cutoff  frequency. 

For  the  examples  to  follow  the  intermediate  covariance  sequence 
approximant  will  be  given,  indicating  the  pole  locations  for  the 
digital  filter.  As  the  denominator  polynomial  A(z)  is  an  intermediate 
result,  it  will  be  given  also.  Determination  of  the  zeros  was  not 
performed  explicitly,  as  different  phase  characteristic  considerations 
will  lead  to  different  polynomials.  Relevant  procedures  for  zero 
selection  are  given  in  Section  5.2. 

Example  1 Desired  stopband  noiselevel:  -80  dB 

In  the  generator  equation  (5.13)  we  use  AN  = 10-^  , corresponding 
to  -40  dB  noiselevel  for  the  intermediate  spectral  density  S^(*)  • 

20  covariance  sequence  lags  were  used; 


Maximum  prod  I c lor  order:  M = 10. 

Lags  7 through  19  determine  the  poles. 

Lags  0 through  8 determine  the  zeros. 

Covariance  Sequence  Approximant : 

Rk  = (.850910"2  - j .86341q'2)(.9774  exp  j .6062)^  + c.c. 

+ (.279010_1  - j .203510_1)(.8932  exp  j .4986)^  + c.c. 

+ (.423910_1  - j . 4219 10'2) (.8087  exp  j .2624) ^ + c.c. 

+ (.416910_1  ) (.8160  )lkl 

A(z)  - 1 - 5.553z_1  + 13. 75z  ' - 19.60z  * + 17.31a"* 

- 9.455z"5  + 2.954z”6  - .4067z"7 


The  resulting  magnitude  squared  characteristic  for  the  recursive 
digital  filter  of  ARMA(14,12)  type  is  given  in  Figure  5.2.  Note  here 
that  the  specifications  practically  call  for  a Pade  like  approximation. 
Cutoff  frequency  peaking  is  one  of  the  characteristics  of  such  an 
approximation.  Pass  band  and  transition  band  performance  are  inferior 
to  stopband  performance  as  might  be  expected  from  the  use  of  short 
covariance  sequences.  Note  also  that  the  maximum  tolerable  system 
order  was  not  used  because  at  lower  system  order  the  prediction  error 
became  negative  indicating  that  the  accuracy  of  the  computer  was  being 
exceeded . 

Example  2 Desired  stopband  noise  level:  -60  dB 

50  covariance  sequence  lags  were  used; 

Maximum  predictor  order:  M - 8. 

Lags  8 through  49  determine  the  poles. 

Lags  0 through  7 determine  the  zeros. 

Covariance  Sequence  Approximant: 


^ - (.283810~2  - j .856010_2)(.9799  exp  j .6223) ^ + c.c. 

+ (.141810_1  + j .233010"1)(.8881  exp  - j .5756) ^ + c.c. 

+ (.443010_1  + j .337310-1)(.7356  exp-1  .4042) ^ + c.c. 

+ (.773810_1  ) (.6401  )lkl 

+ (.778510_3  )(.226710_1  exp  j 

A(z)  - 1.  - 5.053z_1  + 11.56z~2  - 15.32z"3  + 12.62Z-4 

- 6.423z"5  + 1.844z-6  - .2171z“7  - .005948z"8 
The  resulting  magnitude  squared  characteristic  for  the  correspond ing 


ARMA(16 ,14)  digital  filter  is  given  in  Figure  5.3.  The  determination 
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of  the  zeros  is  still  Pade— like  and  accounts  for  the  occurrence  of 
slight  peaking  at  the  cutoff  frequency.  Otherwise  very  reasonable 
performance  is  the  result  in  the  passband,  transition  band  and  stopband. 
Example  3 Desired  stopband  noise  level:  -60  dB 
50  covariance  sequence  lags  were  used; 

Maximum  predictor  order:  M = 10. 


Lags  10  through  49  determine  the  poles. 

Lags  0 through  19  determine  the  zeros. 

Covariance  Sequence  Approximant: 

\ - (.200010-2  - j .6748l0-2)(.9840  exp  j .6239) ^ + c.c. 

+ (.852110_2  + j .180510_1)(.9117  exp-1  .5945)^  + c.c. 

+ (.264610_1  - j .306910_1)(.7811  exp  j .4931) ^ + c.c. 


+ (.627910_1  - j . 2408 1Q-1) (.6427  exp  j .2193) 


+ (-.1209 
+ (.1247 


-5 

10 

-2 


k| 

1 + c.c. 
) (.3537  exp  j ir  ) ^1 


10 


) ( . 2081 


-1 


10 


) 


A(z)  = 1 - 5.406z_1  + 13.22z~2  - 18.67z“3  + 16.11z~4  - 7.987z-5 
+ 1.430z“6  + .6712z'7  - .4401z~8  + .08059z'9 


- .001493z 


-10 


Figure  5.4  indicates  that  a very  nice  result  has  been  achieved  in 
the  passband,  transltionband  and  stopband.  The  use  of  longer  sequences 
has  practically  eliminated  all  ripple  in  the  passband,  but  simultaneously 
the  desired  stopband  noise  level  is  achieved.  The  cutoff  rate  in  the 
transltionband  increases  rapidly  from  about  150  to  600  db/octave. 


Example  4 Desired  stopband  noise  level:  -120  dB 
SO  covariance  sequence  lags  were  used; 

Maximum  predictor  order:  M - 8. 

Lags  8 through  49  determine  the  poles. 

Lags  0 through  19  determine  the  zeros. 

Covariance  Sequence  Approximant : 

R k - (.304210-2  - j .717310_2)(.9829  exp  j .6216)1*1  + c.c. 

+ (.126810_1  - j .176510_1)(.9091  exp  j .5770) '*1  + c.c. 

+ (.321410*1  + j .229910_1)(.7954  exp-)  .4407) '*1  + c.c. 

+ (.521410_1  + j .115510"1)(.7063  exp-j  .1707)^  + c.c. 

A(z)  - 1 - 5.953*"1  + 16.20z"2  - 26.15z-3  + 27.26a"4 

- 18.76z"5  + 8.302*”6  - 2.158*"7  + .2520z-8 
The  magnitude  squared  response  for  the  resulting  ARMA(16,14) 
recursive  digital  filter  is  depicted  in  Figure  5.5.  We  note  that  a 
deterioration  of  the  stopband  performance  takes  place.  Recall  that  the 
stopband  specification  was  merely  an  implicit  design  requirement.  As 
noted  in  Chapter  3 the  underlying  error  criterion  weighs  most  heavily 
the  frequency  regions  with  high  spectral  density,  which  corresponds  to 
the  passband  frequencies.  It  is  therefore  not  surprising  to  see  very 
good  passband  behavior  and  deteriorating  stopband  and  transitionband 
performance.  Even  so  the  initial  fall  off  rate  is  quite  high  and  the 
stopband  rejection  eventually  reaches  -120  dB. 
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Example  5 Desired  stopband  noise  level:  -60  dB 
50  covariance  sequence  lags  were  used; 

Maximum  predictor  order:  M ■ 4. 

Lags  4 through  49  determine  the  poles. 

Lags  0 through  49  determine  the  zeros. 

Covariance  Sequence  Approximant: 

\ =*  (.193110_1  - j .261610_1)(.9383  exp  j .5936) W + c.c. 

+ (.812410_1  + j .194610_1)(.7725  exp-J  .2671) ^ + c.c. 

A(z)  - 1 - 3.064z_1  + 3.796z-2  - 2.240z"3  + .5255s”4 
Figure  5.6  shows  the  magnitude  squared  response  for  this  ABMA(8,6) 
recursive  digital  filter  design.  In  spite  of  the  use  of  relatively 
long  covariance  sequences  the  result  is  not  very  palatable.  The  pass- 
band  has  ripple  of  several  dB,  the  transition  is  not  very  sharp  and  the 
stopband  never  reaches  the  desired  rejection.  This  design  example 
therefore  merely  shows  that  with  too  low  a filter  order  specification 
the  resulting  approximations  cannot  be  expected  to  be  very  good.  The 
latter  of  course  is  true  for  any  filter  design  procedure. 

The  covariance  sequence  approximation  procedure  forms  an  alternative 
method  for  solving  the  stochastic  realization  problem  encountered  in 
the  design  of  covariance  invariant  digital  filters.  The  same  procedure 
may  be  used  to  find  an  approximation  to  an  arbitrary  covariance  sequence 
for  which  no  exact  underlying  AKMA(M,N)  process  exists.  The  designs 
resulting  from  covariance  sequence  approximation  for  modified  ideal 
low  pass  filter  designs  represent  causal,  stable  ARMA(2M,2M~)  recursive 
digital  filters  with  real  coefficients.  We  note  that  absolutely  no 
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phase  Information  is  used  In  the  design.  In  designs  utilizing  Impulse 
response  sequences  at  least  Implicit  use  of  phase  information  is  made. 

A stability  requirement  is  easily  implementable . In  the  designs 
of  this  thesis  an  unstable  result  was  never  encountered,  even  without 
the  use  of  a stability  check.  For  the  Cl  digital  filter  design  the 
covariance  function  for  a stable  continuous  time  process  will  lead  to 
a stable  discrete  time  process  [B] , [PS],  [SP],  but  no  such  proof  is 
available  for  the  approximation  of  arbitrary  covariance  sequences. 
Whereas  no  phase  characteristics  enter  in  the  magnitude  squared  response 
approximation,  considerable  freedom  exists  in  the  numerator  polynomial 
factorization.  As  a result  this  freedom  can  be  used  towards  achieve- 
ment of  linear  phase  characteristics  or  phase  equalization.  Alterna- 
tively minimum  phase  properties  may  be  achieved. 

With  the  covariance  sequence  approximant  procedure  we  therefore 
have  a simple  method  to  design  stable  recursive  filters  of  ARMA(2M,2M  ) 
type  by  solving  linear  systems  of  equations  and  one  polynomial  root- 
finding, both  of  order  M.  The  resulting  designs  compare  favorably  with 
filters  of  the  same  complexity  according  to  a recent  modified  least 
squares  procedure  [MR],  [SL],  even  though  only  half  of  the  filterco- 
efflcients  are  truly  design  variables  due  to  the  doubling  of  poles  and 
zeros  of  the  intermediate  covariance  sequence  approximant.  The  weighted 
least  squares  error  criterion  in  the  modified  least  squares  procedure 
puts  heavy  emphasis  on  the  stopband  performance  and  obviously  at  the 
expense  of  increased  passband  ripple. 

Each  design  was  performed  in  single  precision  on  the  CYBER  171 
at  Colorado  State  University  and  completed  in  less  than  one  second  of 
computation  time.  Due  to  the  fact  that  a spectral  approximation  takes 
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place  one  can  determine  the  magnitude  squared  response  of  the  recursive 
digital  filter  before  actually  determining  its  filter  coefficients.  In 
an  interactive  computation  mode  this  presents  the  possibility  of  a 
visual  check  on  the  resulting  performance  without  actually  performing 
the  complete  design. 
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6 CONCLUSIONS 


6 . 1 Summary 

In  this  thesis  we  have  proposed  a covariance  sequence  model 
consisting  of  a linear  combination  of  damped  complex  exponentials.  A 
generalized  stochastic  process  generator  has  been  derived  that  fits  in 
the  framework  of  the  covariance  sequence  model,  as  its  spectral  repre- 
sentation for  the  covariance  sequence  is  a positive  real  linear  combina- 
tion of  damped  complex  exponentials.  We  have  thus  developed  a generali- 
zation of  the  stochastic  almost  periodic  function  by  allowing  the  complex 
exponentials  to  be  damped.  A more  efficient  representation  for  a large 
subclass  of  ARMA(M,M  ) processes  is  the  result. 

To  identify  the  parameters  in  the  covariance  sequence  model  a 
modified  least  squares  approach  has  resulted  in  a decoupling  of  the 
determination  of  the  complex  exponentials  and  complex  linear  weights. 

Only  linear  systems  of  equations  have  to  be  solved  in  addition  to  a 
polynomial  rootfinding.  A recursive  algorithm  is  given  to  solve  for 
the  complex  exponentials  and  it  is  indicated  how  the  stability  of  the 
resulting  filter  can  be  guaranteed. 

When  used  in  a spectral  analysis  context  the  covariance  sequence 
approximant  of  this  thesis  seems  robust  against  the  nonrat ional 
Gaussian  spectral  density.  For  known  covariance  sequences  with  under- 
lying ARMA(M,M  ) systems,  the  spectral  information  is  obtained  without 
error.  When  a covariance  sequence  is  estimated  from  a finite  data  record 
the  effects  of  the  finite  observation  length  show  up  in  the  covariance 
sequence  estimate  and  consequently  in  the  covariance  sequence  approxi- 
mant.  These  effects  are  derived  and  illustrated  with  realizations  from 
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averaging  do  not  yield  the  same  results  for  the  estimated  covariance 
sequence.  Large  error  terms  can  arise  in  the  identified  complex  linear 
weights.  Error  terms  due  to  difference  frequencies  of  sinusoidal  com- 
ponents can  be  eliminated  by  the  use  of  an  appropriate  length  finite 
observation  record,  resulting  in  accurate  power  and  frequency  estimates 
for  these  signals. 

Harmonic  processes  in  white  noise  form  a limiting  case  for  the 
proposed  covariance  sequence  approximant.  Estimated  covariance  sequences 
for  such  processes  will  yield  relevant  information  concerning  the 
ODserved  process,  but  for  relatively  high  noise  level  the  signal 
components  may  be  averaged  or  not  resolved  at  all. 

In  the  application  of  the  covariance  sequence  approximant  to 
recursive  digital  filter  design  from  magnitude  squared  information  only, 
the  covariance  sequence  approximant  is  sought  for  the  square  root  of  the 
desired  magnitude  squared  response.  This  approach  avoids  by  design  the 
problems  of  an  approximating  spectrum  that  is  possibly  negative  for 
some  frequencies.  The  result  is  a procedure  for  the  design  of  causal 
ARMA(2M,2M  ) filters  with  real  coefficients  of  guaranteed  stability. 

This  design  procedure  requires  line”-  systems  of  equations  to  be  solved 
in  addition  to  a polynomial  rootfinding  problem  and  an  inverse  DFT. 

The  resulting  spectral  error  criterion  for  the  proposed  design  procedure 
is  somewhere  between  modified  least  squares  and  least  squares.  The 
relative  proximity  to  either  is  determined  by  the  choice  of  the  system 
order  and  the  length  of  the  covariance  sequences  used  to  determine 
poles  and  zeros  respectively. 
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6.2  Recommendations 

We  have  seen  that  the  proposed  covariance  sequence  model  is  so 
general  that  approximating  sequences  are  not  necessarily  nonnegative 
definite.  Especially  for  recursive  digital  filter  design  this  is  not 
tolerable.  This  prompted  the  use  of  a square  root  spectrum  approximation 
as  an  intermediate  step.  Finding  general  conditions  for  the  nonnegative' 
ness  of  the  spectral  estimate  could  result  in  a constrained  quadratic 
minimization  problem.  The  solution  then  would  be  an  ARMA(M,M  ) digital 
filter  of  half  the  complexity  of  the  designs  produced  by  Chapter  5.  A 
spectral  factorization  procedure  would  be  unavoidable  however  in  order 
to  find  appropriate  numerator  polynomial  coefficients. 

The  covariance  sequence  approxlmant  of  this  thesis  is  valid  under 
the  assumption  of  single  poles  for  the  underlying  system.  A different 
covariance  sequence  approxlmant  can  be  defined  to  incorporate  multiple 
poles.  In  the  solution  procedure  this  would  result  in  derivatives  of 
the  independent  vectors  in  the  Vandermonde  matrix.  One  can  argue  that 
practical  systems  identified  under  noisy  conditions  will  hardly  ever 
exhibit  multiple  poles.  The  importance  of  this  generalized  approach 
then  lies  in  the  fact  that  multiple  poles  at  zero  can  be  implemented 
which  will  allow  the  identification  of  ARMA(M,N)  systems,  where  N 
can  be  larger  than  M-l. 

While  illustrative,  the  spectrum  analysis  results  of  Chapter  4 are 
not  conclusive.  Continued  investigation  should  focus  on  Monte-Carlo 
studies  of  spectrum  estimates  for  random  amplitude,  random  phase 
sinusoids,  together  with  analytical  investigation  of  estimator  variances, 
confidence  intervals,  and  the  like. 
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